##%%%%%%%%%%%%%% simOfDeming functions %%%%%%%%%%%%%%%%%%%%%%
## To get the number of occurances of a value in a 
## a = table(ordinalVar)
## We will use this to compute the lists of distributions and supports
## needed for the functions in the GenOrd package

calcSupportM = function(XX){
	## This function assumes XX is a matrix of ordinal variables
	## it then computes lists of marginal distributions and supports
	## for these using the table function as well as the correlation matrix.
	supportList=list()
	marginalList=list()
	varNames=colnames(XX)
	for(i in 1:ncol(XX)){
		tabI=table(XX[,i])
		supportList[[i]]=as.numeric(names(tabI))
		margI=tabI/sum(tabI)
		margI=cumsum(margI)[1:(length(tabI)-1)]
		marginalList[[i]]=margI	
		}
	names(supportList)=varNames
	names(marginalList)=varNames
	corXX=cor(XX)
	outList=list()
	outList[[1]]=marginalList
	outList[[2]]=corXX
	outList[[3]]=supportList
	names(outList)=c("marginal", "Sigma", "support")
	outList
	}


## It may turn out that the correlation matrix produced here won't be
## feasible.  We then use the following function to bring it into range.

supportCorrFix=function(simInput){
	## Use 'corrcheck' function from GenOrd package to get
	## upper and lowerbounds for correlation matrix and replace
	## values in the correlation matrix with feasible values.
 	corBounds=corrcheck(marginal=simInput$marginal,support=simInput$support)
	# Lower bounds
	corL=corBounds[[1]]
	# Upper bounds
	corU=corBounds[[2]]
	corM=simInput$Sigma
	violL=as.matrix(corM<corL)
	violU=as.matrix(corM>corU)
	corOut=corM
	corOut[violL]=corL[violL]
	corOut[violU]=corU[violU]
	simInputOut=simInput
	simInputOut$Sigma=corOut
 	simInputOut
	}

## Here we use the histogram function in r to establish binsides
## and counts

makeOrdinal = function(x){
	## Here if there are less than 21 unique
	## elements in the vector we leave it alone
	## otherwise we use the histogram function
	## to compute bins, probabilities
	ouT=x
	tab1=table(x)
	if(length(tab1)>20){
		hh=hist(x)
		N=length(hh$mids)
		breakS=hh$breaks
		midS=hh$mids
 		for(i in 1:(N-1)){
			breakSi=breakS[i:(i+1)]
			midSi=midS[i]
	 		iinn= (x>=breakSi[1]) & (x<breakSi[2])
	 		ouT[iinn]=midSi
			}
		iiEnd=(x>=breakS[N])
		ouT[iiEnd]=midS[N]
		}
	ouT
	}

## Matrix version of this function
makeOrdinalM = function(XX){
 	## XX should be a matrix on which
	## we do column wise "makeOrdinal"
	ouT=XX
 	for(i in 1:ncol(XX)){
		ouT[,i]=makeOrdinal(XX[,i])
		}
 	ouT
	}




## Multi-binning function with plotting

Unif = function(x,ncats=9){
	JJ = is.finite(x)
	X=x[JJ]
	binsides = unique(quantile(X,(0:ncats)/ncats))
	if(length(binsides)<3){stop("one bin or less in Unif")}
	condMat = matrix(T,length(X),ncats)
	condMat[,1] = (X >= binsides[1] & X <= binsides[2])
	for(i in 3:length(binsides)){ condMat[,i-1] = (X > binsides[i-1] & X <= binsides[i])}
	ouT = x
	valS = as(2*((1:ncats) -mean(1:ncats))/(ncats-1),"matrix")
	ouT[JJ] = condMat%*%valS
	ouT
}

MultiBin = function(x,y,ncats = 6,
		binSides=c(),
		yNames=c("Y"),
        xName="X",
		rounD=4,
		rcor="unif",
		maiN=""){
	## this function bins the x variable, then computes 
	## mean values and stdev of each of the y-variables
	## conditioned on the x variable being in bin-i
	## Function assumes length(x) matches nrows-y
	## Check to make sure that dimensions are OK
	if(length(x) != nrow(as(y,"matrix"))){ stop("dim x and y don't match")}
	X = junk2zero(x)
	Y = junk2zero(as(y,"matrix"))
    if(length(yNames)<(dim(Y)[2])){ yNames=paste("Y",1:((dim(Y)[2])),sep="")}
	## compute binsides for x
	if(length(binSides) > 0){ binsides = binSides}
	if(length(binSides) == 0){ binsides = unique(quantile(X,(0:ncats)/ncats))}
	if(length(binsides)<3){stop("one bin or less")}
	
	## compute conditioning matrix
	condMat = matrix(0,length(X),length(binsides)-1)
	condMat[,1] = (X >= binsides[1] & X <= binsides[2])
	for(i in 3:length(binsides)){ condMat[,i-1] = (X > binsides[i-1] & X <= binsides[i])}
	
	## test for empty bins
	foo = Csum(condMat)
	condMat = condMat[,foo > 0]
	binsides = c(binsides[1],binsides[2:length(binsides)][foo>0])
	
	## Compute conditional means and stdev's
	norM = junk2zero(1/t(matrix(1,(dim(Y)[2]),1)%*%Csum(condMat)))
	YmeanS = t(t(Y)%*%condMat)*norM
	YstdevS = sqrt(t(t(Y^2)%*%condMat)*norM -YmeanS^2)*sqrt(norM)
		
	## Compute x-means and stdev's
	norM = junk2zero(1/Csum(condMat))
	XmeanS = (t(as(X,"matrix"))%*%condMat)*norM
	XstdevS = sqrt((t(as(X^2,"matrix"))%*%condMat)*norM - XmeanS^2)*sqrt(norM)

	## produce labels 
	xLab = paste(round(binsides[1:(length(binsides)-1)],rounD),"<x<=",
	                          round(binsides[2:length(binsides)],rounD),sep="")
	dimnames(YmeanS) = list(xLab,paste("cond-mean ",yNames,sep=""))
	dimnames(YstdevS) = list(xLab,paste("cond std ",yNames,sep=""))
	
	## compute correlations
	corS = cor(X,Y)
	rcorS = corS-corS
	if(rcor=="unif"){for(i in 1:(dim(Y)[2])){rcorS[i] = cor(Unif(X,10),Unif(Y[,i],10))}}
	if(rcor=="rank"){for(i in 1:(dim(Y)[2])){rcorS[i] = cor(rank(X),rank(Y[,i]))}}

	corS = as(round(corS,3),"numeric")
	rcorS = as(round(rcorS,3),"numeric")
	
	names(corS) = yNames
	names(rcorS) = yNames
	
	## output
	
	ouT = list(yMeans = YmeanS,yStdevs = YstdevS,xMeans = XmeanS,
	                    xStdevs=XstdevS,cors=corS,rcors=rcorS,
						yNames=yNames,xName=xName,maiN=maiN)
	ouT
}
	

## MultiBin using bin medians
MultiBinM = function(x,y,ncats = 6,
		binSides=c(),
		yNames=c("Y"),
        xName="X",
		rounD=4,
		rcor="unif",
		maiN=""){
	## this function bins the x variable, then computes 
	## mean values and stdev of each of the y-variables
	## conditioned on the x variable being in bin-i
	## Function assumes length(x) matches nrows-y
	## Check to make sure that dimensions are OK
	if(length(x) != nrow(as(y,"matrix"))){ stop("dim x and y don't match")}
	X = junk2zero(x)
	Y = junk2zero(as(y,"matrix"))
    if(length(yNames)<(dim(Y)[2])){ yNames=paste("Y",1:((dim(Y)[2])),sep="")}
	## compute binsides for x
	if(length(binSides) > 0){ binsides = binSides}
	if(length(binSides) == 0){ binsides = unique(quantile(X,(0:ncats)/ncats))}
	if(length(binsides)<3){stop("one bin or less")}
	
	## compute conditioning matrix
	condMat = matrix(0,length(X),length(binsides)-1)
	condMat[,1] = (X >= binsides[1] & X <= binsides[2])
	for(i in 3:length(binsides)){ condMat[,i-1] = (X > binsides[i-1] & X <= binsides[i])}
	
	## test for empty bins
	foo = Csum(condMat)
	condMat = condMat[,foo > 0]
	binsides = c(binsides[1],binsides[2:length(binsides)][foo>0])
	
	## Compute conditional means and stdev's
	norM = junk2zero(1/t(matrix(1,(dim(Y)[2]),1)%*%Csum(condMat)))
	YmeanS = t(t(Y)%*%condMat)*norM
	for(i in 1:length(YmeanS)){YmeanS[i,1]=median(Y[as.logical(condMat[,i])])}
	YstdevS = sqrt(t(t(Y^2)%*%condMat)*norM -YmeanS^2)*sqrt(norM)
		
	## Compute x-means and stdev's
	norM = junk2zero(1/Csum(condMat))
	XmeanS = (t(as(X,"matrix"))%*%condMat)*norM
	XstdevS = sqrt((t(as(X^2,"matrix"))%*%condMat)*norM - XmeanS^2)*sqrt(norM)

	## produce labels 
	xLab = paste(round(binsides[1:(length(binsides)-1)],rounD),"<x<=",
	                          round(binsides[2:length(binsides)],rounD),sep="")
	dimnames(YmeanS) = list(xLab,paste("cond-median ",yNames,sep=""))
	dimnames(YstdevS) = list(xLab,paste("cond-std ",yNames,sep=""))
	
	## compute correlations
	corS = cor(X,Y)
	rcorS = corS-corS
	if(rcor=="unif"){for(i in 1:(dim(Y)[2])){rcorS[i] = cor(Unif(X,10),Unif(Y[,i],10))}}
	if(rcor=="rank"){for(i in 1:(dim(Y)[2])){rcorS[i] = cor(rank(X),rank(Y[,i]))}}

	corS = as(round(corS,3),"numeric")
	rcorS = as(round(rcorS,3),"numeric")
	
	names(corS) = yNames
	names(rcorS) = yNames
	
	## output
	
	ouT = list(yMeans = YmeanS,yStdevs = YstdevS,xMeans = XmeanS,
	                    xStdevs=XstdevS,cors=corS,rcors=rcorS,
						yNames=yNames,xName=xName,maiN=maiN)
	ouT
}
	
	
binPlot = function(xMean,yMean,yStdev=sqrt(var(yMean)),stdS=2,yName="Y",xName="X",cor=0,rcor=0,maiN=""){
		yMat = matrix(0,3,length(yMean))
		yMat[2,] = yMean
		yMat[1,] = yMean-yStdev*stdS
		yMat[3,] = yMean+yStdev*stdS
		rangE = max(yMat)-min(yMat)
		yliM = c(min(yMat)-0.05*rangE,max(yMat)+0.05*rangE)
        Main=maiN
		if(maiN=="") Main=paste("condMean ",yName," v. ",xName,",",stdS,"std-err")
        ##par(cex=0.5)
		plot(xMean,yMean,ylim=yliM,main=Main,
		      xlab=paste("x=",xName,"      cor=",cor," rcor=",rcor,sep=""),ylab=yName)
        ## par(cex=1)
		for(i in 1:length(xMean)){
			lines(c(1,1,1)*xMean[i],yMat[,i])
		}
      }



CerrPlot=function(y,xMean=(0:(ncol(y)-1)),maiN="",yName="Y",xName="X",lineS=F,stdErr=T){
  yMat=matrix(0,3,ncol(y))
  yMean=Cmean(y)
  if(stdErr==T) yStdev=2*Cstd(y)/sqrt(Csum(sign(is.finite(y))))
  if(stdErr==F) yStdev=2*Cstd(y)
  yMat[2,] = yMean
  yMat[1,] = yMean-yStdev
  yMat[3,] = yMean+yStdev
  rangE = max(yMat)-min(yMat)
  yliM = c(min(yMat)-0.05*rangE,max(yMat)+0.05*rangE)
  Main=maiN
  if(maiN=="") Main=paste("Mean ",yName," v. ",xName,", 2-std-err")
  if(lineS==F){plot(xMean,yMean,ylim=yliM,main=Main,ylab=yName,xlab=xName)
               for(i in 1:length(xMean)){
                 lines(c(1,1,1)*xMean[i],yMat[,i])
               }
             }
  if(lineS==T){yliM[1]=min(yliM[1],0)
               yliM[2]=max(yliM[2],0)
               plot(c(0-0.5,xMean),c(0,yMean),ylim=yliM,main=Main,ylab=yName,xlab=xName)
               for(i in 1:length(xMean)){
                 lines(c(1,1,1)*xMean[i],yMat[,i])
               }
               lines(c(0-0.5,xMean),c(0,yMean),type="b")
             }
  
}


RerrPlot=function(y,xMean=(0:(ncol(t(y))-1)),maiN="",yName="Y",xName="Y"){
  yMat=matrix(0,3,ncol(t(y)))
  yMean=Cmean(t(y))
  yStdev=2*Cstd(t(y))/sqrt(Csum(sign(is.finite(t(y)))))
  yMat[2,] = yMean
  yMat[1,] = yMean-yStdev
  yMat[3,] = yMean+yStdev
  rangE = max(yMat)-min(yMat)
  yliM = c(min(yMat)-0.05*rangE,max(yMat)+0.05*rangE)
  Main=maiN
  if(maiN=="") Main=paste("Mean ",yName," v. ",xName,", 2-std-err")
  plot(xMean,yMean,ylim=yliM,main=Main,ylab=yName,xlab=xName)
  for(i in 1:length(xMean)){
    lines(c(1,1,1)*xMean[i],yMat[,i])
  }
}


MultiBinPlot = function(mb,stdS=2,mfrowRevert=T,lineS=F){
	## This is a plotting function for the output of MultiBin
	if(ncol(mb$yMeans)>1){
	nplotS = ceiling(ncol(mb$yMeans)/2)
	par(mfrow=c(nplotS,2))
}
	for(i in 1:ncol(mb$yMeans)){
		binPlot(mb$xMeans,mb$yMeans[,i],mb$yStdevs[,i],
		        stdS=stdS,yName=mb$yNames[i],xName=mb$xName,cor=mb$cors[i],rcor=mb$rcors,maiN=mb$maiN)
	}
	if(mfrowRevert==T & ncol(mb$yMeans)>1) par(mfrow=c(1,1))
    if(lineS==T) lines(mb$xMeans,mb$xMeans)
}
	

	
ErrPlot = function(x,y,ncats= 6,
		binSides=c(),
		yName=c("Y"),
        xName=c("X"),
		rounD=4,
		rcor="unif",
		maiN="",
        lineS=F){
			Y = as(y,"matrix")
			foobar = MultiBin(x,Y,ncats=ncats,binSides=binSides,yNames=yName,xName=xName,rounD=rounD,rcor=rcor,maiN=maiN)
			MultiBinPlot(foobar,lineS=lineS)
			foobar
		}
	

CondPlot = function(x,y,criT,ncats=6,critCats=6,
							binSides=c(),
							yName=c("Y"),
							rounD=4,
							rcor="unif",
							maiN="",
							lineS=F){
				Y = as(y,"matrix")
				qq=unique(quantile(criT,(0:critCats)/critCats))
				nPlots=length(qq)-1
				if(nPlots < 2) stop("Only one condition catagory!")
				par(mfrow=c(ceiling(nPlots/2),2))
				ouT=list()
                for(i in 1:(length(qq)-1)){
					cc=criT>qq[i] & criT<=qq[i+1]
					foobar = MultiBin(x[cc],Y[cc,],ncats=ncats,binSides=binSides,yNames=yName,rounD=rounD,rcor=rcor,maiN=paste(maiN,round(qq[i],rounD),"<z<",round(qq[i+1],rounD),sep=""))
					MultiBinPlot(foobar,mfrowRevert=F)
					if(lineS==T) lines(foobar$xMeans,foobar$xMeans)
					ouT[[i]]=foobar
					}
		par(mfrow=c(1,1))
		ouT
       }




eventPlot = function(yMean,yStdev,timeS=c(0,1:length(yMean)),yName="y",stdS=2,xbin="xbin"){
	## Helper plotting function for event plotting
	if(length(timeS)!=length(yMean)+1 || sum(timeS==0)!=1){stop("timeS wrong in eventPlot")}
	zeroIndex = (1:length(yMean))[timeS==0]
	yMat = matrix(0,3,length(timeS))
	ii = 2:length(timeS)
	if(zeroIndex > 1 && zeroIndex <length(timeS)){
		ii = c(1:(zeroIndex-1),(zeroIndex+1):length(timeS))
		}
	if(zeroIndex==length(timeS)){ii = 1:(length(timeS)-1)}
	yMat[1,ii] = yMean-yStdev*stdS
	yMat[2,ii] = yMean
	yMat[3,ii] = yMean+yStdev*stdS
	rangE = max(yMat)-min(yMat)
	yliM = c(min(yMat)-0.025*rangE,max(yMat)+0.025*rangE)
	plot(timeS,yMat[2,],main=paste(yName,"over time","for",xbin),type="l",ylim=yliM,ylab=yName)
	lines(timeS,yMat[1,], lty=8)
	lines(timeS,yMat[3,],lty=8)
	lines(timeS,yMat[1,]-yMat[1,])
}


MultiBinEventPlot = function(mb,timeS=c(0,1:ncol(mb$yMeans)),stdS=2,yName="y"){
	## This is a plotting function for event analysis. It assumes a MultiBin object 
	## has been created that is conditioned around an event. The dependent variables 
	## of the MultiBin object are assumed to be of returns or something similar. Times
	## refer to times around the event. The zero in times gives the index where the
	## look-back times of the event change to look forward.
	
	## Check that timeS has the right length and has the other correct properties
	if(length(timeS) != ncol(mb$yMeans)+1){stop("timeS wrong length!")}
	if(sum(timeS[2:length(timeS)] > timeS[1:(length(timeS)-1)]) < length(timeS)-1){
		stop("timeS not monotonly increasing")}
	if(sum(timeS == 0)==0){stop("timeS doesn't contain a zero!")}
	nplotS = ceiling(length(mb$xMeans)/2)
	par(mfrow=c(nplotS,2))
	xbinS = dimnames(mb$yMeans)[[1]]
	for(i in 1:length(mb$xMeans)){
		eventPlot(yMean=mb$yMeans[i,],yStdev=mb$yStdevs[i,],timeS=timeS,
		                yName=yName,stdS=stdS,xbin=xbinS[i])
	}
	par(mfrow=c(1,1))
	
}


	
## Testing

lagM = function(x,i=0){
	ouT = as(x,"matrix")
	if(i>0){
		ouT[(i+1):nrow(ouT),] = ouT[1:(nrow(ouT)-i),]
		ouT[1:i,] = 0
	}
	if(i < 0){
		ouT[1:(nrow(ouT)-abs(i)),] = ouT[(abs(i)+1):nrow(ouT),]
		ouT[(nrow(ouT)-abs(i)+1):nrow(ouT),] = 0
	}
	ouT
}

winSum = function(x,i){
	ouT = as(x,"matrix")
	if(i > 0){
		for(l in 1:i) ouT = ouT+lagM(x,l)
	}
	if(i < 0){
		ouT = ouT-ouT
		for(l in -1:i) ouT = ouT+lagM(x,l)
	}
	ouT	
}

corErr = function(coR,nVal){ouT = (1-coR^2)/sqrt(nVal)
           ouT}

bootCor = function(x, y, n = 100)
{
	bootMat = matrix(0, n, 1)
	for(i in 1:n) {
		ii = sample(1:length(x),replace=T)
		bootMat[i] = cor(x[ii], y[ii])
	}
	qN = c(0.9,0.95,0.975)
	ouT = round(quantile(bootMat, qN),3)
	ouT
}

## boot correlation for matrices
bootCorM = function(X,Y,n=100){
  bootMat=matrix(0,n,ncol(X)*ncol(Y))
  for(i in 1:n) {
		ii = sample(1:nrow(X),replace=T)
		bootMat[i,] = as.vector(cor(X[ii,], Y[ii,]))
	}
  errVec=as.vector(Cstd(bootMat))
  errMat=matrix(errVec,ncol(X),ncol(Y))
  dimnames(errMat)=list(dimnames(X)[[2]],dimnames(Y)[[2]])
  errMat}


resSub = function(x,II =list(1:ncol(x))){
	ouT = x
	for(i in 1:length(II)){
		ii = II[[i]]
		ouT[,ii] = ouT[,ii] - (Rmean(x[,ii])%*%matrix(1,1,ncol(x[,ii])))
			}
	ouT
}






##%%%%%%%% Row Column operators  %%%%%%%%%%%%%%%%%%%

Rsum = function(x){
	ouT = as(junk2zero(x)%*%matrix(1,ncol(x),1),"matrix")
	dimnames(ouT) = list(dimnames(x)[[1]],c("sum"))
	ouT
}


Csum = function(x){
	ouT = as(matrix(1,1,nrow(x))%*%junk2zero(x),"matrix")
	dimnames(ouT) = list(c("sum"),dimnames(x)[[2]])
	ouT
}

Rmean = function(x){
	ouT = as(junk2zero(x)%*%matrix(1,ncol(x),1),"matrix")*junk2zero(1/(is.finite(x)%*%matrix(1,ncol(x),1)))
	dimnames(ouT) = list(dimnames(x)[[1]],c("sum"))
	ouT
}


Cmean = function(x){
	ouT = (as(matrix(1,1,nrow(x))%*%junk2zero(x),"matrix")*(junk2zero(1/(matrix(1,1,nrow(x))%*%is.finite(x)))))
	dimnames(ouT) = list(c("sum"),dimnames(x)[[2]])
	ouT
}

CmeanNZ = function(x){
	ouT = (as(matrix(1,1,nrow(x))%*%junk2zero(x),"matrix")*(junk2zero(1/(matrix(1,1,nrow(x))%*%(is.finite(x)&(x!=0))))))
	dimnames(ouT) = list(c("sum"),dimnames(x)[[2]])
	ouT
}

Cstd = function(x){
	meanS = matrix(1,nrow(x),1)%*%Cmean(x)
	ouT = junk2zero(sqrt(Cmean((x-meanS)^2)))
	dimnames(ouT) = list(c("sum"),dimnames(x)[[2]])
	ouT
}

Rstd = function(x){
	meanS = Rmean(x)%*%matrix(1,1,ncol(x))
	ouT = junk2zero(sqrt(Rmean((x-meanS)^2)))
	dimnames(ouT) = list(dimnames(x)[[1]],c("sum"))
	ouT
}

csumMat = function(j){
	II = (1:j)%*%matrix(1,1,j)
	JJ = matrix(1,j,1)%*%(1:j)
	ouT = (JJ >= II)+0
}

CcumSum = function(x){
	ouT = x
	if(ncol(x) > 100){stop("number of cols greater than 100!")}	
	for(i in 1:ncol(x)){ouT[,i] = cumsum(x[,i])}
	ouT
}


RcumSum = function(x){
	ouT = junk2zero(x)
	if(ncol(x) > 500){stop("number of cols greater than 500!")}
   ouT = ouT%*%csumMat(ncol(ouT))
	ouT
}

Rmax = function(x){
	ouT=junk2zero(x[,1])
	if(ncol(x) > 150) {stop("number of cols greater than 150!")}
	for(i in 2:ncol(x)){ ouT[junk2zero(x[,i])>ouT]=x[junk2zero(x[,i])>ouT,i]}
	ouT}

Rmin = function(x){
	ouT=junk2zero(x[,1])
	if(ncol(x) > 150) {stop("number of cols greater than 150!")}
	for(i in 2:ncol(x)){ ouT[junk2zero(x[,i])<ouT]=x[junk2zero(x[,i])<ouT,i]}
	ouT}

Cmax = function(x){
	if(ncol(x) > 150) {stop("number of cols greater than 150!")}
	ouT=junk2zero(x[1,])
	for(i in 1:ncol(x)){ ouT[i] = max(junk2zero(x[,i]))}
	ouT}


Cmin = function(x){
	if(ncol(x) > 150) {stop("number of cols greater than 150!")}
	ouT=junk2zero(x[1,])
	for(i in 1:ncol(x)){ ouT[i] = min(x[is.finite(x[,i]),i])}
	ouT}


## With NA handling
Rmaxs = function(x){
	if(ncol(x) > 500) {stop("number of cols greater than 500!")}
	z=x
 	ii=!is.finite(z)
	naVal=min(z[!ii])-1
  	z[ii]=naVal
	ouT=z[,1]
	for(i in 2:ncol(x)){ ouT[z[,i]>ouT]=z[z[,i]>ouT,i]}
	jj=ouT==naVal
	ouT[jj]=NA
	ouT}

Rmins = function(x){
	if(ncol(x) > 500) {stop("number of cols greater than 500!")}
	z=x
 	ii=!is.finite(z)
	naVal=max(z[!ii])+1
  	z[ii]=naVal
	ouT=z[,1]
	for(i in 2:ncol(x)){ ouT[z[,i]<ouT]=z[z[,i]<ouT,i]}
	jj=ouT==naVal
	ouT[jj]=NA
	ouT}

Cmaxs = function(x){
	if(ncol(x) > 500) {stop("number of cols greater than 500!")}
	z=x
 	ii=!is.finite(z)
	naVal=min(z[!ii])-1
  	z[ii]=naVal
      ouT=z[1,]
	for(i in 1:ncol(x)){ ouT[i] = max(z[,i])}
	jj=ouT==naVal
	ouT[jj]=NA
	ouT}


Cmins = function(x){
	if(ncol(x) > 500) {stop("number of cols greater than 500!")}
	z=x
 	ii=!is.finite(z)
	naVal=max(z[!ii])+1
  	z[ii]=naVal
	ouT=z[1,]
	for(i in 1:ncol(x)){ ouT[i] = min(z[,i])}
	jj=ouT==naVal
	ouT[jj]=NA
	ouT}



Rmeans=function(M){
 	s1=Rsum(j2z(M))
	n1=Rsum(!is.na(M))
 	ii=n1>0
 	ouT=j2z(M[,1]-M[,1])
 	ouT[ii]=s1[ii]/n1[ii]
	ouT[!ii]=NA
	ouT}

Rsums=function(M){
 	s1=Rsum(j2z(M))
	n1=Rsum(!is.na(M))
 	ii=n1>0
 	ouT=j2z(M[,1]-M[,1])
 	ouT[ii]=s1[ii]
	ouT[!ii]=NA
	ouT}

Csums=function(M){
      ouT=t(Rsums(t(M)))
      ouT}

Cmeans=function(M){
      ouT=t(Rmeans(t(M)))
      ouT}

Rstds=function(M){
      meanM=Rmeans(M)%*%matrix(1,1,ncol(M))
      dM=(M-meanM)^2
      dMM=Rsums(dM)
      norM=Rsums(!is.na(M))-1
      norM[norM==0]=1
      dMMM=dMM/norM
      ouT=sqrt(dMMM)
      ouT}

Cstds=function(M){
      ouT=t(Rmeans(t(M)))
      ouT}




stds=function(x){
 ii=!is.na(x)
 mm=mean(x[ii])
 ss=sqrt(var(x[ii]))
 ouT=(x-mm)/ss
 ouT}

Cmedian = function(x){
	if(ncol(x) > 150) {stop("number of cols greater than 150!")}
	ouT=junk2zero(x[1,])
	for(i in 1:ncol(x)){ ouT[i] = median(x[is.finite(x[,i]),i])}
	ouT}

CmedianNZ = function(x){
	if(ncol(x) > 150) {stop("number of cols greater than 150!")}
	ouT=junk2zero(x[1,])
	for(i in 1:ncol(x)){ ouT[i] = median(x[is.finite(x[,i])&(x[,i]!=0),i])}
	ouT}

##%%%%%%%%%% pair functions %%%%%%%%%%%%%%%

 Ccor=function(x, y)
	{
	xMean = Cmean(x)
	yMean = Cmean(y)
	xSigma = sqrt(Cmean((x - matrix(1,nrow(y),1)%*%xMean)^2))
	ySigma = sqrt(Cmean((y - matrix(1,nrow(y),1)%*%yMean)^2))
	xyMean = Cmean(x * y)
	rcor = (xyMean - xMean * yMean)/(xSigma * ySigma)
	rcor
	}

 Rcor= function(x, y)
	{
	xMean = Rmean(x)
	yMean = Rmean(y)
	xSigma = sqrt(Rmean((x - xMean %*% matrix(1, 1, ncol(x)))^2))
	ySigma = sqrt(Rmean((y - yMean %*% matrix(1, 1, ncol(x)))^2))
	xyMean = Rmean(x * y)
	rcor = (xyMean - xMean * yMean)/(xSigma * ySigma)
	rcor
	}



##%%%%%%%% Edge window stuff ######################


RedgeWin=function(x,edgE){
  NN=unique(Rsum(sign(edgE)))
  if(length(NN)>1){ stop("window edges inconsistent!")}
  ## Flip things to make it work
  ouT=t(matrix(t(x)[t(edgE)],NN,nrow(x)))
  ouT}

CedgeWin=function(x,edgE){
  NN=unique(Csum(sign(edgE)))
  if(length(NN)>1){ stop("window edges inconsistent!")}
  ouT=matrix(x[edgE],NN,ncol(x))
  ouT}
  
##%%%%%%%%%%% Array functions for 3d array on last dim %%%%%%%%%%

Zsum = function(x){
    zzn=dim(x)[[3]]
	if(zzn > 150) {stop("number of dim-3 greater than 100!")}
	ouT=junk2zero(x[,,1])
	for(i in 2:zzn){ ouT = ouT+junk2zero(x[,,i])}
	ouT}

Zmean = function(x){
    norM=junk2zero((1/Zsum(sign(is.finite(x)))))
	ouT = Zsum(x)*norM
	ouT
}

Zmin = function(x){
    zzn=dim(x)[[3]]
	if(zzn > 150) {stop("number of dim-3 greater than 100!")}
	ouT=junk2zero(x[,,1])
	for(i in 2:zzn){ ouT = miN(ouT,junk2zero(x[,,i]))}
	ouT}

Zmax = function(x){
    zzn=dim(x)[[3]]
	if(zzn > 150) {stop("number of dim-3 greater than 100!")}
	ouT=junk2zero(x[,,1])
	for(i in 2:zzn){ ouT = maX(ouT,junk2zero(x[,,i]))}
	ouT}

Zstd = function(x){
    zzn=dim(x)[[3]]
	if(zzn > 150) {stop("number of dim-3 greater than 100!")}
	meanS = x-x
    foo=Zmean(x)
    for(i in 1:zzn){meanS[,,i]=foo}
	ouT = junk2zero(sqrt(Zmean((junk2zero(x)-meanS)^2)))
	ouT}

minN0=function(x,y){
  ouT=x
  ouT[(y<x)&(y>0)]=y[(y<x)&(y>0)]
  ouT}

Zmin0 = function(x){
    zzn=dim(x)[[3]]
	if(zzn > 150) {stop("number of dim-3 greater than 100!")}
	ouT=junk2zero(x[,,1])
	for(i in 2:zzn){ ouT = miN0(ouT,junk2zero(x[,,i]))}
	ouT}



##%%%%%%%%%% Functions to do aggregations "by" an index vector %%%%%%%

sumBy = function(x,idx){
  idxU=unique(idx)
  ouT=matrix(0,length(idxU),1)
  names(ouT)=idxU
  for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i]=sum(x[jj])}
  ouT}


meanBy = function(x,idx){
  idxU=unique(idx)
  ouT=matrix(0,length(idxU),1)
  names(ouT)=idxU
  for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i]=mean(x[jj])}
  ouT}

medianBy = function(x,idx){
  idxU=unique(idx)
  ouT=matrix(0,length(idxU),1)
  names(ouT)=idxU
  for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i]=median(x[jj])}
  ouT}



minBy = function(x,idx){
  idxU=unique(idx)
  ouT=matrix(0,length(idxU),1)
  names(ouT)=idxU
  for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i]=min(x[jj])}
  ouT}

maxBy = function(x,idx){
  idxU=unique(idx)
  ouT=matrix(0,length(idxU),1)
  names(ouT)=idxU
  for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i]=max(x[jj])}
  ouT}


stdBy = function(x,idx){
  idxU=unique(idx)
  ouT=matrix(0,length(idxU),1)
  names(ouT)=idxU
  for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i]=sqrt(var(x[jj]))}
  ouT}

firstBy=function(x, idx){
	idxU = unique(idx)
	ouT = matrix(0, length(idxU), 1)
	names(ouT) = idxU
	for(i in 1:length(idxU)) {
		jj = (idx == idxU[i])
		ouT[i] = x[jj][1]
	}
	ouT}


lastBy=function(x, idx){
	idxU = unique(idx)
	ouT = matrix(0, length(idxU), 1)
	names(ouT) = idxU
	for(i in 1:length(idxU)) {
		jj = (idx == idxU[i])
		ouT[i] = x[jj][sum(jj)]
	}
	ouT}

## The same row and column-wise

RsumBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[,1:length(idxU)]
    dimnames(ouT)=list(dimnames(x)[[1]],idxU)
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[,i]=Rsum(x[,jj])}
    ouT}
      

CsumBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[1:length(idxU),]
    dimnames(ouT)=list(idxU,dimnames(x)[[2]])
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i,]=Csum(x[jj,])}
    ouT}


RmeanBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[,1:length(idxU)]
    dimnames(ouT)=list(dimnames(x)[[1]],idxU)
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[,i]=Rmean(x[,jj])}
    ouT}
      

CmeanBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[1:length(idxU),]
    dimnames(ouT)=list(idxU,dimnames(x)[[2]])
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i,]=Cmean(x[jj,])}
    ouT}

CmeanByNZ = function(x,idx){
    idxU=unique(idx)
    ouT=x[1:length(idxU),]
    dimnames(ouT)=list(idxU,dimnames(x)[[2]])
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i,]=CmeanNZ(x[jj,])}
    ouT}


RstdBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[,1:length(idxU)]
    dimnames(ouT)=list(dimnames(x)[[1]],idxU)
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[,i]=Rstd(x[,jj])}
    ouT}
      

CstdBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[1:length(idxU),]
    dimnames(ouT)=list(idxU,dimnames(x)[[2]])
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i,]=Cstd(x[jj,])}
    ouT}


RmaxBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[,1:length(idxU)]
    dimnames(ouT)=list(dimnames(x)[[1]],idxU)
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[,i]=Rmax(x[,jj])}
    ouT}
      

CmaxBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[1:length(idxU),]
    dimnames(ouT)=list(idxU,dimnames(x)[[2]])
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i,]=Cmax(x[jj,])}
    ouT}

RminBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[,1:length(idxU)]
    dimnames(ouT)=list(dimnames(x)[[1]],idxU)
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[,i]=Rmin(x[,jj])}
    ouT}
      

CminBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[1:length(idxU),]
    dimnames(ouT)=list(idxU,dimnames(x)[[2]])
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i,]=Cmin(x[jj,])}
    ouT}


CmedianBy = function(x,idx){
    idxU=unique(idx)
    ouT=x[1:length(idxU),]
    dimnames(ouT)=list(idxU,dimnames(x)[[2]])
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i,]=Cmedian(x[jj,])}
    ouT}


CmedianByNZ = function(x,idx){
    idxU=unique(idx)
    ouT=x[1:length(idxU),]
    dimnames(ouT)=list(idxU,dimnames(x)[[2]])
    for(i in 1:length(idxU))
      {jj=(idx==idxU[i])
       ouT[i,]=CmedianNZ(x[jj,])}
    ouT}



##%%%%%%%%%%%%%%%%%%%% Trim, NA replace stuff %%%%%%%%%%%%%%%%%
## Trim Tails

TrimW = function(x,triM=0.0,replacE = T){
  trM = quantile(abs(x),1-triM/2)
  CC = is.finite(x) & abs(x) >= trM
	ouT = x
	if(replacE == T){ouT[CC] = sign(ouT[CC])*trM} else {ouT[CC] = 0}
		ouT		
}


## Trim inside quantile values. Sets inside values (abs(x)<qq) to zero
## and pulls outside values so that they go out from zero. (where abs(x)>qq,
## newx=oldx-qq*sign(oldx)
TrimI = function(x,triM=0.0){
	qq=quantile(abs(x[is.finite(x)]),triM)
	II = is.finite(x) & abs(x) <=qq
	TT = is.finite(x) & abs(x) > qq
	ouT = x
	ouT[II] = 0
	ouT[TT]=ouT[TT]-qq*sign(ouT[TT])
		ouT		
}



junk2zero=function(x)
  {
        ouT = x
        if(class(x) == "data.frame") {
                for(i in 1:ncol(ouT)) {
                        if(is(ouT[1, i], "numeric") == T) {
                                ouT[!is.finite(ouT[, i]), i] = 0
                        }
                }
        }
        else ouT[!is.finite(ouT)] = 0
        ouT
   }

j2z=junk2zero

Unif = function(x,ncats=9){
	JJ = is.finite(x)
	X=x[JJ]
	binsides = unique(quantile(X,(0:ncats)/ncats))
	if(length(binsides)<3){stop("one bin or less in Unif")}
	condMat = matrix(T,length(X),ncats)
	condMat[,1] = (X >= binsides[1] & X <= binsides[2])
	for(i in 3:length(binsides)){ condMat[,i-1] = (X > binsides[i-1] & X <= binsides[i])}
	ouT = x
	valS = as(2*((1:ncats) -mean(1:ncats))/(ncats-1),"matrix")
	ouT[JJ] = condMat%*%valS
	ouT
}


qqMat= function(x,ncats=9){
	JJ = is.finite(x)
	X=x[JJ]
	binsides = unique(quantile(X,(0:ncats)/ncats))
	if(length(binsides)<3){stop("one bin or less in Unif")}
	condMat = matrix(T,length(X),ncats)
	condMat[,1] = (X >= binsides[1] & X <= binsides[2])
	for(i in 3:length(binsides)){ condMat[,i-1] = (X > binsides[i-1] & X <= binsides[i])}
	ouT = matrix(F,length(x),ncats)
	ouT[JJ,] = condMat
	ouT
}


fillS=function(x){
  ouT=x
  firstX=((1:length(x))[is.finite(x)])[1]
  ouT[1:firstMp]=x[firstMp]
  jj=cumsum(sign(is.finite(ouT)))
  outStart=ouT[is.finite(ouT)]
  ouTT=outStart[jj]
  ouTT}


##%%%%%%%%%%%%% Discretize continuous variables %%%%%%%%%%%%%%%%%%%
## Here we use the histogram function in r to establish binsides
## and counts

makeOrdinal = function(x){
	## Here if there are less than 21 unique
	## elements in the vector we leave it alone
	## otherwise we use the histogram function
	## to compute bins, probabilities
	ouT=x
	tab1=table(x)
	if(length(tab1)>20){
		hh=hist(x)
		N=length(hh$mids)
		breakS=hh$breaks
		midS=hh$mids
 		for(i in 1:(N-1)){
			breakSi=breakS[i:(i+1)]
			midSi=midS[i]
	 		iinn= (x>=breakSi[1]) & (x<breakSi[2])
	 		ouT[iinn]=midSi
			}
		iiEnd=(x>=breakS[N])
		ouT[iiEnd]=midS[N]
		}
	ouT
	}

## Matrix version of this function
makeOrdinalM = function(XX){
 	## XX should be a matrix on which
	## we do column wise "makeOrdinal"
	ouT=XX
 	for(i in 1:ncol(XX)){
		ouT[,i]=makeOrdinal(XX[,i])
		}
 	ouT
	}


## calcSupportM is used to calculate the various elements
## needed for the functions in the GenOrd package
## To get the number of occurances of a value in a 
## a = table(ordinalVar)
## We will use this to compute the lists of distributions and supports
## needed for the functions in the GenOrd package

calcSupportM = function(XX){
	## This function assumes XX is a matrix of ordinal variables
	## it then computes lists of marginal distributions and supports
	## for these using the table function as well as the correlation matrix.
	supportList=list()
	marginalList=list()
	varNames=colnames(XX)
	for(i in 1:ncol(XX)){
		tabI=table(XX[,i])
		supportList[[i]]=as.numeric(names(tabI))
		margI=tabI/sum(tabI)
		margI=cumsum(margI)[1:(length(tabI)-1)]
		marginalList[[i]]=margI	
		}
	names(supportList)=varNames
	names(marginalList)=varNames
	corXX=cor(XX)
	outList=list()
	outList[[1]]=marginalList
	outList[[2]]=corXX
	outList[[3]]=supportList
	names(outList)=c("marginal", "Sigma", "support")
	outList
	}


## It may turn out that the correlation matrix produced here won't be
## feasible.  We then use the following function to bring it into range.

supportCorrFix=function(simInput){
	## Use 'corrcheck' function from GenOrd package to get
	## upper and lowerbounds for correlation matrix and replace
	## values in the correlation matrix with feasible values.
 	corBounds=corrcheck(marginal=simInput$marginal,support=simInput$support)
	# Lower bounds
	corL=corBounds[[1]]
	# Upper bounds
	corU=corBounds[[2]]
	corM=simInput$Sigma
	violL=as.matrix(corM<corL)
	violU=as.matrix(corM>corU)
	corOut=corM
	corOut[violL]=corL[violL]
	corOut[violU]=corU[violU]
	simInputOut=simInput
	simInputOut$Sigma=corOut
 	simInputOut
	}


##%%%%%%%%%%%%%%%%%%%% Cohen's f^2 %%%%%%%%%%%%%%%%%%%%%%%%%%
## Define function to compute Cohen's f^2 for a set of variables

fsquared = function(tar,explanVars,indexVec){
	nameS=colnames(explanVars)[indexVec]
	fsqAll=matrix(0,length(indexVec),1)
	rownames(fsqAll)=nameS
	colnames(fsqAll)=c("f^2")
 	RsqAll=summary(lm(tar ~ explanVars))$r.squared
 	for(i in 1:length(indexVec)){
		Idrop=setdiff(1:ncol(explanVars),indexVec[i])
		explanVdrop=explanVars[,Idrop]
		RsqDrop=summary(lm(tar ~ explanVars[,Idrop]))$r.squared
		fsqAll[i,]=(RsqAll-RsqDrop)/(1-RsqAll)
		}
	fsqAll
	}


##%%%%%%%%%% Function to convert Fstat to tstat %%%%%%%%%%%%%

FtoT=function(Fstat,df1,df2){
	ouT=Fstat-Fstat
	jj=Fstat>0
	## Go from 1-side F-prob to 
	## two sided t-prob
	Fprob=pf(Fstat[jj],df1,df2)
	F2Tprob=1-(1-Fprob)/2
	ouT[jj]=qt(F2Tprob,df2)
	ll=!is.finite(ouT[jj])
	ouT[jj][ll]=4.6+0.0541*Fstat[jj][ll]
	ouT
	}


##%%%%%%%%%%%%% Odds and Ends %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
norM=function(vec){ouT=(vec-mean(vec))/sqrt(var(vec)); ouT}

rsqOut=function(taR,inS,coeF){
	coef=coeF
	coef[!is.finite(coef)]=0
	Vartar=var(taR)
	fitVal=coef[1]+inS%*%(as.matrix(coef[2:length(coef)]))
	Varres=var(taR-fitVal)
	rsQ=1-Varres/Vartar
	rsQ
	}


adj.rsqOut=function(taR,inS,coeF){
	coef=coeF
	n=length(taR)
	coef[!is.finite(coef)]=0
	p=sum(coef!=0)-1
	adj=(n-1)/(n-p-1)
	Vartar=var(taR)
	fitVal=coef[1]+inS%*%(as.matrix(coef[2:length(coef)]))
	Varres=var(taR-fitVal)
	rsQ=1-(Varres*adj)/Vartar
	rsQ
	}


Rmedian=function(M){
	ouT=matrix(0,nrow(M),1)
	rownames(ouT)=rownames(M)
	colnames(ouT)=c("row median")
	for(i in 1:nrow(M)){
		ouT[i,]=median(M[i,])
		}
	ouT 	
	}

Cmedian=function(M){
	ouT=matrix(0,1,ncol(M))
	colnames(ouT)=colnames(M)
	rownames(ouT)=c("column median")
	for(i in 1:ncol(M)){
		ouT[,i]=median(M[,i])
		}
	ouT 	
	}


Rquantile=function(M,qq){
	ouT=matrix(0,nrow(M),length(qq))
	rownames(ouT)=rownames(M)
	for(i in 1:nrow(M)){
		rqq=quantile(M[i,],qq)
		ouT[i,]=rqq
		}
	colnames(ouT)=names(rqq)
	ouT
	}

Cquantile=function(M,qq){
	ouT=matrix(0,length(qq),ncol(M))
	colnames(ouT)=colnames(M)
	for(i in 1:ncol(M)){
		cqq=quantile(M[,i],qq)
		ouT[,i]=cqq
		}
	rownames(ouT)=names(cqq)
	ouT
	}



##%%%%%%%%%% Compute in-sample model lists %%%%%%%%%%%%
## Function to compute modList
makeModList=function(tar,inS,IndexMat){
	N=ncol(IndexMat)
	modList=list()
	for(i in 1:N){
		ii=IndexMat[,i]
		t1=tar[ii]
		ins=inS[ii,]
		modI=lm(t1 ~ ins)
		modIsum=summary(modI)
		shortList=list()
		shortList[[1]]=modI
		shortList[[2]]=modIsum
		names(shortList)=c("modI","modIsum")
		modList[[i]]=shortList
		}
	names(modList)=paste("mod",1:N,sep="")
	modList
	}

## Function to extract matrix of model coefficients
## from a modList

makeCoefMat=function(modList){
	N=length(modList)
	coef1=modList[[1]][[1]]$coefficients
	coefmat=matrix(0,length(coef1),N)
	rownames(coefmat)=names(coef1)
	colnames(coefmat)=names(modList)
	for(i in 1:N){
		coefmat[,i]=modList[[i]][[1]]$coefficients
		}
	coefmat[!is.finite(coefmat)]=0
	coefmat
	}

## Function to extract matrix of coefficient tstats
## from a modList

makeTstatMat=function(modList){
	N=length(modList)
	coef1=modList[[1]][[1]]$coefficients
	coefmat=matrix(0,length(coef1),N)
	rownames(coefmat)=names(coef1)
	colnames(coefmat)=names(modList)
	for(i in 1:N){
		coefSum=modList[[i]][[2]]$coefficients
		coefNames=rownames(coefSum)
		coefmat[coefNames,i]=coefSum[coefNames,3]
		}
	coefmat[!is.finite(coefmat)]=0
	coefmat
	}
	
## Function to extract array of coefficient stats
## from a modList

makeCoefArr=function(modList){
	N=length(modList)
	coef1=modList[[1]][[1]]$coefficients
	coefSum1=modList[[1]][[1]]$coefficients
	coefarr=array(0,dim=c(length(coef1),N,4))
	dimnames(coefarr)=list(names(coef1),names(modList),colnames(coefSum1))
	for(i in 1:N){
		coefSum=modList[[i]][[2]]$coefficients
		coefNames=rownames(coefSum)
		coefarr[coefNames,i,]=coefSum[coefNames,]
		}
	coefarr[!is.finite(coefarr)]=0
	coefarr
	}


##%%%%%%%% Make rsq list function with cross validation %%%%%%%%%%%
## For this function one needs to call "library(MASS)"
## for the Moore-Penrose psuedo inverse "ginv"

crossValD1=function(taR,Ins,lmModel,sumModel){
	iiS=is.finite(lmModel$coef[2:length(lmModel$coef)])
	X=cbind(matrix(1,nrow(Ins),1),Ins)
	hatMat=X%*%ginv(t(X)%*%X)%*%t(X)
	trHatMat=diag(hatMat)
	trHatMat[trHatMat>0.5]=0.5
	resiD=sumModel$residuals
	norM=1/(1-trHatMat)
	norM[!is.finite(norM)]=0
	## The following correction is to 
	## prevent numerical outliers where
	## the trHatMat~1 making norM's 
	## very large.
	norM[norM>2]=2
	cvRes=resiD*norM
	vaRcvRes=var(cvRes)

	rsqMat=matrix(0,4,1)
	rownames(rsqMat)=c("r.squared","adj.r.squared",
					"cv.r.squared","cv.adj.r.squared")
	n=length(taR)
	p=sum(iiS)
	adj=(n-1)/(n-p-1)
	rsqMat[1,]=sumModel$r.squared
	rsqMat[2,]=sumModel$adj.r.squared
	rsqMat[3,]=1-vaRcvRes/var(taR)
	rsqMat[4,]=1-(adj*vaRcvRes)/var(taR)
	rsqMat
	}



## This version of the function doesn't require
## lmModel or sumModel to have been pre-computed

crossValD1S=function(taR,Ins){
	lmModel=lm(taR~Ins)
	iiS=is.finite(lmModel$coef[2:length(lmModel$coef)])
	sumModel=summary(lmModel)
	X=cbind(matrix(1,nrow(Ins),1),Ins)
	hatMat=X%*%ginv(t(X)%*%X)%*%t(X)
	trHatMat=diag(hatMat)
	trHatMat[trHatMat>0.5]=0.5
	resiD=sumModel$residuals
	norM=1/(1-trHatMat)
	norM[!is.finite(norM)]=0
	norM[norM>2]=2
	cvRes=resiD*norM
	vaRcvRes=var(cvRes)

	rsqMat=matrix(0,4,1)
	rownames(rsqMat)=c("r.squared","adj.r.squared",
					"cv.r.squared","cv.adj.r.squared")
	n=length(taR)
	p=sum(iiS)
	adj=(n-1)/(n-p-1)
	rsqMat[1,]=sumModel$r.squared
	rsqMat[2,]=sumModel$adj.r.squared
	rsqMat[3,]=1-vaRcvRes/var(taR)
	rsqMat[4,]=1-(adj*vaRcvRes)/var(taR)
	rsqMat
	}

## This is the traditional leave one out CV

crossValB1=function(tar,inS){
	lmModel=lm(tar~inS)
	iiS=is.finite(lmModel$coef[2:length(lmModel$coef)])
	inSS=inS[,iiS]
	cvyHat=matrix(0,ncol(inSS),1)
	for(i in 1:nrow(inSS)){
		jjd1=setdiff(1:nrow(inSS),i)
		cv.coef=j2z(lm(tar[jjd1]~inSS[jjd1,])$coef)
		cvyHat[i]=cv.coef[1]+inSS[i,]%*%as.matrix(cv.coef[2:length(cv.coef)])
		}
	cvRes=tar-cvyHat
	cv.r.squared=1-var(cvRes)/var(tar)
	cv.r.squared
	}





##%%%%%%%% CrossVal function for a set of inputs %%%%%%%%%%%
crossValD=function(taR,Ins,jjVec){
	cv.rsqMat=matrix(0,6,length(jjVec))
	colnames(cv.rsqMat)=colnames(Ins[,jjVec])				
	rownames(cv.rsqMat)=c("r.squaredI","r.squaredD","cv.r.squaredI",
					"cv.r.squaredD","f.squared","cv.f.squared")

	modI=lm(taR~Ins)
	modIsum=summary(modI)
	cvRes=crossValD1(taR,Ins,modI,modIsum)
	cv.rsqMat[1,]=cvRes[1]
	cv.rsqMat[3,]=cvRes[3]
	
	for(i in 1:length(jjVec)){
		jjD=setdiff(1:ncol(Ins),jjVec[i])
		modD=lm(taR~Ins[,jjD])
		modDsum=summary(modD)
		cvResD=crossValD1(taR,Ins[,jjD],modD,modDsum)
		cv.rsqMat[2,i]=cvResD[1]
		cv.rsqMat[4,i]=cvResD[3]
		}
	cv.rsqMat[5,]=(cv.rsqMat[1,]-cv.rsqMat[2,])/(1-cv.rsqMat[1,])
	cv.rsqMat[6,]=(cv.rsqMat[3,]-cv.rsqMat[4,])/(1-cv.rsqMat[3,])
	cv.rsqMat
	}


## Traditional leave one out cross validation for a set of inputs
crossValB=function(taR,Ins,jjVec){
	cv.rsqMat=matrix(0,6,length(jjVec))
	colnames(cv.rsqMat)=colnames(Ins[,jjVec])				
	rownames(cv.rsqMat)=c("r.squaredI","r.squaredD","cv1.r.squaredI",
					"cv1.r.squaredD","f.squared","cv1.f.squared")

	modI=lm(taR~Ins)
	modIsum=summary(modI)
	cvRes=crossValB1(taR,Ins,modI,modIsum)
	cv.rsqMat[1,]=cvRes[1]
	cv.rsqMat[3,]=cvRes[3]
	
	for(i in 1:length(jjVec)){
		jjD=setdiff(1:ncol(Ins),jjVec[i])
		modD=lm(taR~Ins[,jjD])
		modDsum=summary(modD)
		cvResD=crossValD1(taR,Ins[,jjD],modD,modDsum)
		cv.rsqMat[2,i]=cvResD[1]
		cv.rsqMat[4,i]=cvResD[3]
		}
	cv.rsqMat[5,]=(cv.rsqMat[1,]-cv.rsqMat[2,])/(1-cv.rsqMat[1,])
	cv.rsqMat[6,]=(cv.rsqMat[3,]-cv.rsqMat[4,])/(1-cv.rsqMat[3,])
	cv.rsqMat
	}




##%%%%%% M-fold cross validation %%%%

## M-fold cross validation for a group of inputs
crossValM=function(tar,inputS,jjVec,M){
	N=nrow(inputS)
	Nout=floor(N/M)	
	iiSm=1:Nout
	partM=(0:(M-1))*Nout
	outList=list()

	for(i in 1:(M-1)){outList[[i]]=iiSm+partM[i]}
	outList[[M]]=(partM[M]+1):N
	cv.rsqMat=matrix(0,6,length(jjVec))
	rownames(cv.rsqMat)=c("r.squaredI","r.squaredD","cv.r.squaredI",
					"cv.r.squaredD","f.squared","cv.f.squared")
	colnames(cv.rsqMat)=colnames(inputS)[jjVec]
	resVecI=tar-tar
	resMatD=matrix(0,length(resVecI),length(jjVec))
	modIsum=summary(lm(tar~inputS))
	cv.rsqMat[1,]=modIsum$r.squared

	for(i in 1:M){
		iiOut=outList[[i]]
		iiIn=setdiff(1:N,iiOut)
		modIcoef=j2z(lm(tar[iiIn]~inputS[iiIn,])$coef)
		resVecI[iiOut]=modIcoef[1]+inputS[iiOut,]%*%as.matrix(modIcoef[2:length(modIcoef)])-tar[iiOut]
		}

	for(j in 1:length(jjVec)){
		jjD=setdiff(1:ncol(inputS),jjVec[j])
		modDsum=summary(lm(tar~inputS[,jjD]))
		cv.rsqMat[2,j]=modDsum$r.squared
		for(i in 1:M){
			iiOut=outList[[i]]	
			iiIn=setdiff(1:N,iiOut)
			modDcoef=j2z(lm(tar[iiIn]~inputS[iiIn,jjD])$coef)
			resMatD[iiOut,j]=modDcoef[1]+inputS[iiOut,jjD]%*%as.matrix(modDcoef[2:length(modDcoef)])-tar[iiOut]
			}
		}
	cv.rsqMat[3,]=1-var(resVecI)/var(tar)
	cv.rsqMat[4,]=1-as.vector(diag(var(resMatD)))/as.numeric(var(tar))
	cv.rsqMat[5,]=(cv.rsqMat[1,]-cv.rsqMat[2,])/(1-cv.rsqMat[1,])
	cv.rsqMat[6,]=(cv.rsqMat[3,]-cv.rsqMat[4,])/(1-cv.rsqMat[3,])
	cv.rsqMat
	}
		


## Function to compute M-fold cross validation for an 
## entire set of models.  This function is designed to have 
## an output that can be easily merged into the results 
## of the "makeRsqList" type functions.

makeCvmList=function(tar,inputS,ImaT,jjVec,M){
	Nmod=ncol(ImaT)
	cvmArray=array(0,dim=c(6,length(jjVec),Nmod))
	for(l in 1:Nmod){
		iiIn=ImaT[,l]
		t1In=tar[iiIn]
		insIn=inputS[iiIn,]
		cvM=crossValM(t1In,insIn,jjVec,M)
		cvmArray[,,l]=cvM
		}
	rsqList=list()
	cvrNames=rownames(cvM)
	cvcNames=paste("mod",1:Nmod,sep="")
		
	for(j in 1:length(jjVec)){
		rsqMat=matrix(0,4,Nmod)
		fsqMat=matrix(0,2,Nmod)
		rownames(rsqMat)=cvrNames[1:4]
		rownames(fsqMat)=cvrNames[5:6]
		colnames(rsqMat)=cvcNames
		colnames(fsqMat)=cvcNames

		for(l in 1:Nmod){
			outVec=cvmArray[,j,l]			
			rsqMat[,l]=outVec[1:4]
			fsqMat[,l]=outVec[5:6]
			}
		shortList=list()
		shortList[[1]]=rsqMat
		shortList[[2]]=fsqMat
		names(shortList)=c("rsqMat","fsqMat")
		rsqList[[j]]=shortList
		}
	names(rsqList)=colnames(inputS)[jjVec]
	rsqList
	}


## Function to compute M-fold cross validation for an 
## entire set of models.  This function is designed to have 
## an output that can be easily merged into the results 
## of the "makeRsqList" type functions.

makeCvmList=function(tar,inputS,ImaT,jjVec,M){
	Nmod=ncol(ImaT)
	cvmArray=array(0,dim=c(6,length(jjVec),Nmod))
	for(l in 1:Nmod){
		iiIn=ImaT[,l]
		t1In=tar[iiIn]
		insIn=inputS[iiIn,]
		cvM=crossValM(t1In,insIn,jjVec,M)
		cvmArray[,,l]=cvM
		}
	rsqList=list()
	cvrNames=rownames(cvM)
	cvcNames=paste("mod",1:Nmod,sep="")
		
	for(j in 1:length(jjVec)){
		rsqMat=matrix(0,4,Nmod)
		fsqMat=matrix(0,2,Nmod)
		rownames(rsqMat)=cvrNames[1:4]
		rownames(fsqMat)=cvrNames[5:6]
		colnames(rsqMat)=cvcNames
		colnames(fsqMat)=cvcNames

		for(l in 1:Nmod){
			outVec=cvmArray[,j,l]			
			rsqMat[,l]=outVec[1:4]
			fsqMat[,l]=outVec[5:6]
			}
		shortList=list()
		shortList[[1]]=rsqMat
		shortList[[2]]=fsqMat
		names(shortList)=c("rsqMat","fsqMat")
		rsqList[[j]]=shortList
		}
	names(rsqList)=colnames(inputS)[jjVec]
	rsqList
	}


## Function to compute simple version of leave one out
## cross validation using trace of hat matrix for an 
## entire set of models.  This function is designed to have 
## an output that can be easily merged into the results 
## of the "makeRsqList" type functions.

makeCvList=function(tar,inputS,ImaT,jjVec){
	Nmod=ncol(ImaT)
	cvmArray=array(0,dim=c(6,length(jjVec),Nmod))
	for(l in 1:Nmod){
		iiIn=ImaT[,l]
		t1In=tar[iiIn]
		insIn=inputS[iiIn,]
		cvM=crossValD(t1In,insIn,jjVec)
		cvmArray[,,l]=cvM
		}
	rsqList=list()
	cvrNames=rownames(cvM)
	cvcNames=paste("mod",1:Nmod,sep="")
		
	for(j in 1:length(jjVec)){
		rsqMat=matrix(0,4,Nmod)
		fsqMat=matrix(0,2,Nmod)
		rownames(rsqMat)=cvrNames[1:4]
		rownames(fsqMat)=cvrNames[5:6]
		colnames(rsqMat)=cvcNames
		colnames(fsqMat)=cvcNames

		for(l in 1:Nmod){
			outVec=cvmArray[,j,l]			
			rsqMat[,l]=outVec[1:4]
			fsqMat[,l]=outVec[5:6]
			}
		shortList=list()
		shortList[[1]]=rsqMat
		shortList[[2]]=fsqMat
		names(shortList)=c("rsqMat","fsqMat")
		rsqList[[j]]=shortList
		}
	names(rsqList)=colnames(inputS)[jjVec]
	rsqList
	}


##%%%%%%%%% Function to compute bootstrap smooth all cv, cvM %%%%%


makeCvaBoot=function(tar,inputS,ImaT,bootMat,jjVec){
	Nmod=ncol(ImaT)
	Nboot=ncol(bootMat)
	cvArray=array(0,dim=c(6,length(jjVec),Nmod,Nboot))
	dimnames(cvArray)[[3]]=paste("mod",1:Nmod,sep="")
	dimnames(cvArray)[[4]]=paste("boot",1:Nboot,sep="")
	dimnames(cvArray)[[2]]=names(jjVec)
	cv10Array=cvArray
	cv5Array=cvArray
	cv2Array=cvArray
	dimnames(cvArray)[[1]]= c("r.squaredI","r.squaredD","cv.r.squaredI",
						"cv.r.squaredD","f.squared","cv.f.squared")
	dimnames(cv10Array)[[1]]= c("r.squaredI","r.squaredD","cv10.r.squaredI",
						"cv10.r.squaredD","f.squared","cv10.f.squared")
	dimnames(cv5Array)[[1]]= c("r.squaredI","r.squaredD","cv5.r.squaredI",
						"cv5.r.squaredD","f.squared","cv5.f.squared")
	dimnames(cv2Array)[[1]]= c("r.squaredI","r.squaredD","cv2.r.squaredI",
						"cv2.r.squaredD","f.squared","cv2.f.squared")
	
	for(j in 1:Nboot){
		iib=bootMat[,j]
		
		for(l in 1:Nmod){
			iiIn=ImaT[iib,l]
			t1In=tar[iiIn]
			insIn=inputS[iiIn,]
			cv=crossValD(t1In,insIn,jjVec)
			cv10=crossValM(t1In,insIn,jjVec,10)
			cv5=crossValM(t1In,insIn,jjVec,5)
			cv2=crossValM(t1In,insIn,jjVec,2)
			cvArray[,,l,j]=cv
			cv10Array[,,l,j]=cv10
			cv5Array[,,l,j]=cv5
			cv2Array[,,l,j]=cv2
			}
		}
	arrayList=list()
	arrayList[[1]]=cvArray
	arrayList[[2]]=cv10Array
	arrayList[[3]]=cv5Array
	arrayList[[4]]=cv2Array
	names(arrayList)=c("cvArray","cv10Array","cv5Array","cv2Array")
	arrayList
	}




##%%%%%%%%%%%%%%% Create merged lists %%%%%%%%%%%%%%%%%%%%%%%%%
## The following function merges original rsqList
## with cross validation lists

mergeLists=function(rsqList,cvL,cvL10,cvL5,cvL2){
	cvList=list()
	for(i in 1:length(cvL10)){
		cv.rsqMatI=rbind(cvL[[i]][[1]][3,],cvL10[[i]][[1]][3,],cvL5[[i]][[1]][3,],
							cvL2[[i]][[1]][3,])
		rownames(cv.rsqMatI)=paste("cv",c("","10","5","2"),".r.squared",sep="")
		cv.rsqMatD=rbind(cvL[[i]][[1]][4,],cvL10[[i]][[1]][4,],cvL5[[i]][[1]][4,],
							cvL2[[i]][[1]][4,])
		rownames(cv.rsqMatD)=paste("cv",c("","10","5","2"),".r.squaredD",sep="")
		cv.fsqMat=rbind(cvL[[i]][[2]][2,],cvL10[[i]][[2]][2,],cvL5[[i]][[2]][2,],
							cvL2[[i]][[2]][2,])
		rownames(cv.fsqMat)=paste("cv",c("","10","5","2"),".f.squared",sep="")
		shortList=list()
		shortList[[1]]=cv.rsqMatI
		shortList[[2]]=cv.rsqMatD
		shortList[[3]]=cv.fsqMat
		cvList[[i]]=shortList
		}

	## Merge the cvLists with rsqList
	mergeList=rsqList
	for(i in 1:length(mergeList)){
		rsqImat=rbind(rsqList[[i]][[1]][c(1:2,5),],cvList[[i]][[1]])
		mergeList[[i]][[1]]=rsqImat
		rsqDmat=rbind(rsqList[[i]][[2]][c(1:2,5),],cvList[[i]][[2]])
		mergeList[[i]][[2]]=rsqDmat
		fsqMat=rbind(rsqList[[i]][[3]][c(1:2,5),],cvList[[i]][[3]])
		mergeList[[i]][[3]]=fsqMat
		}
	mergeList
	}
	



## The following function replaces the old cv.data
## with the properly computed cv.data

mergeListCV=function(cv.rsqList,cvL1){
	cvList=cv.rsqList
	for(i in 1:length(cvL1)){
		cvList[[i]][[1]]["cv.r.squared",]=cvL1[[i]][[1]]["cv.r.squaredI",]
		cvList[[i]][[2]]["cv.r.squared",]=cvL1[[i]][[1]]["cv.r.squaredD",]
		cvList[[i]][[3]]["cv.f.squared",]=cvL1[[i]][[2]]["cv.f.squared",]
		}
	cvList
	}
	


##%%%%%%% Function to Merge cvA.rsq lists with bootCvaM lists %%%%%%%%%%%

mergeRsqLBoot=function(cvA.rsq,outL.rsq,outL.fsq,bootCvaM){
	cvl.rsq=cvA.rsq
	for(i in 1:length(cvA.rsq)){
		varN=names(cvA.rsq)[i]
		rsqI.cvA=cvA.rsq[[varN]][[1]]
		rsqD.cvA=cvA.rsq[[varN]][[2]]
		fsq.cvA=cvA.rsq[[varN]][[3]]

		rsqI.outL=outL.rsq[1,]
		rsqD.outL=outL.rsq[varN,]
		fsq.outL=outL.fsq[varN,]

		rsqI.bootM=rbind(bootCvaM[[1]][3,varN,],bootCvaM[[2]][3,varN,],
					bootCvaM[[3]][3,varN,],bootCvaM[[2]][3,varN,])
		rsqD.bootM=rbind(bootCvaM[[1]][4,varN,],bootCvaM[[2]][4,varN,],
					bootCvaM[[3]][4,varN,],bootCvaM[[2]][4,varN,])
		fsq.bootM=rbind(bootCvaM[[1]][6,varN,],bootCvaM[[2]][6,varN,],
					bootCvaM[[3]][6,varN,],bootCvaM[[2]][6,varN,])

		rsqIall=rbind(rsqI.cvA[1:3,],rsqI.outL,rsqI.cvA[4:7,],rsqI.bootM)
		rsqDall=rbind(rsqD.cvA[1:3,],rsqD.outL,rsqD.cvA[4:7,],rsqD.bootM)
		fsqall=rbind(fsq.cvA[1:3,],fsq.outL,fsq.cvA[4:7,],fsq.bootM)
		
		preFix=c("","adj.","out.","outL.","cv.","cv10.","cv5.","cv2.",
				"cvb.","cv10b.","cv5b.","cv2b.")
		rownames(rsqIall)=paste(preFix,"r.squaredI",sep="")
		rownames(rsqDall)=paste(preFix,"r.squaredD",sep="")
		rownames(fsqall)=paste(preFix,"f.squared",sep="")

		cvl.rsq[[i]][[1]]=rsqIall
		cvl.rsq[[i]][[2]]=rsqDall
		cvl.rsq[[i]][[3]]=fsqall
		}
	cvl.rsq
	}



## Function to merge just the long out of sample
mergeRsqL=function(cvA.rsq,outL.rsq,outL.fsq){
	cvl.rsq=cvA.rsq
	for(i in 1:length(cvA.rsq)){
		varN=names(cvA.rsq)[i]
		rsqI.cvA=cvA.rsq[[varN]][[1]]
		rsqD.cvA=cvA.rsq[[varN]][[2]]
		fsq.cvA=cvA.rsq[[varN]][[3]]

		rsqI.outL=outL.rsq[1,]
		rsqD.outL=outL.rsq[varN,]
		fsq.outL=outL.fsq[varN,]

		rsqIall=rbind(rsqI.cvA[1:3,],rsqI.outL,rsqI.cvA[4:7,])
		rsqDall=rbind(rsqD.cvA[1:3,],rsqD.outL,rsqD.cvA[4:7,])
		fsqall=rbind(fsq.cvA[1:3,],fsq.outL,fsq.cvA[4:7,])
		
		preFix=c("","adj.","out.","outL.","cv.","cv10.","cv5.","cv2.")
		rownames(rsqIall)=paste(preFix,"r.squaredI",sep="")
		rownames(rsqDall)=paste(preFix,"r.squaredD",sep="")
		rownames(fsqall)=paste(preFix,"f.squared",sep="")

		cvl.rsq[[i]][[1]]=rsqIall
		cvl.rsq[[i]][[2]]=rsqDall
		cvl.rsq[[i]][[3]]=fsqall
		}
	cvl.rsq
	}





##%%%%%% Function to compute r^2, f^2, anova-stats etc. %%%%%%%

makeRsqList=function(tar,inS,Imat,ImatOut,jjVec){
	jjList=(1:length(jjVec))
	names(jjList)=as.character(jjVec)
	rsqList=list()

    for(j in jjVec){
	Nmod=ncol(Imat)
	jjj=1:ncol(inS)
	jj=setdiff(jjj,j)
	jjjcoef=1:(ncol(inS)+1)
	jjcoef=setdiff(jjjcoef,j+1)
	coefMatD=matrix(0,length(jjcoef),Nmod)
	rownames(coefMatD)=c("intercept",colnames(inS)[jj])
	colnames(coefMatD)=paste("mod",1:Nmod,sep="")
	rsqVecI=matrix(0,4,Nmod)
	rownames(rsqVecI)=c("r.squared","adj.r.squared","out.r.squared","out.adj.r.squared")
	rsqVecD=rsqVecI
	fsqIdealVec=matrix(0,1,Nmod)
	rownames(fsqIdealVec)=c("fsqIdeal")

	## residual coefficients are computed on the residuals
	## of the dropped models
	resCoefVec=matrix(0,2,Nmod)
	rownames(resCoefVec)=c("coef","resCoefInD")
	

	## 	n=number of obs; p=number of inputs;
	##	df_modI=(n-p-1); df modD=(n-(p-1)-1)=(n-p)
	##	F-stat  = ((RSS_modD - RSS_modI)/(df_modD=df_modI)) / (RSS_modI/df_modI)
	##          = ((var(res_modD)-var(res_modI))*(n-1)) / ((var(res_modI)*(n-1))/(n-p-1))
	##    Pr(>F)  = 1 - pf(F-stat,df_modD1-df_modS,df_modS)
	
	anovaMatIn=matrix(0,6,Nmod)
	rownames(anovaMatIn)=c("df_modI","RSS_modD","RSS_modI","SumOfSqD","F","Pr(>F)")
	anovaMatOut=anovaMatIn

	for(i in 1:Nmod){
		iiin=Imat[,i]
		t1in=tar[iiin]
		insI=inS[iiin,]
		insD=inS[iiin,jj]
		modD=lm(t1in ~ insD)
		modI=lm(t1in ~insI)
		modDsum=summary(modD)
		modIsum=summary(modI)
		coefD=modD$coefficients
		coefD[!is.finite(coefD)]=0
		coefI=modI$coefficients
		coefI[!is.finite(coefI)]=0
		coefMatD[,i]=coefD
		rsqVecD[1,i]=modDsum$r.squared
		rsqVecD[2,i]=modDsum$adj.r.squared
		rsqVecI[1,i]=modIsum$r.squared
		rsqVecI[2,i]=modIsum$adj.r.squared
		
		
		## Compute residual coefficient
		resD=modD$residuals
		inD=inS[iiin,j]
		modRes=lm(resD~inD)
		resCoefVec[2,i]=modRes$coefficients[2]
		resCoefVec[1,i]=j2z(coefI[j+1])

		## Compute the out of sample r.squared, r.squared.drop
		iiout=ImatOut[,i] ## Note ImatOut
		t1out=tar[iiout]
		insDout=inS[iiout,jj]
		insIout=inS[iiout,]
		rsqVecD[3,i]=rsqOut(t1out,insDout,coefD)
		rsqVecD[4,i]=adj.rsqOut(t1out,insDout,coefD)
		rsqVecI[3,i]=rsqOut(t1out,insIout,coefI)
		rsqVecI[4,i]=adj.rsqOut(t1out,insIout,coefI)

		
		## Compute fsqIdeal
		v1=inS[iiin,j]
		coef1=coefI[j+1]
		fsqI=(var(v1*coef1)/var(t1in))/(1-rsqVecI[1,i])
		fsqIdealVec[i]=fsqI
	
		## Compute anova SS and F values
		anovaD=anova(modD,modI)
		n=length(iiout);p=sum(j2z(coefI)!=0)
		 # In-sample c("df_modI","RSS_modD","RSS_modI","SumOfSqD","F","Pr(>F)")
		anovaMatIn["df_modI",i]=anovaD[["Res.Df"]][2]
		anovaMatIn[c("RSS_modD","RSS_modI"),i]=anovaD[["RSS"]]
		anovaMatIn["SumOfSqD",i]=anovaD[["Sum of Sq"]][2]
		anovaMatIn["F",i]=anovaD[["F"]][2]
		anovaMatIn["Pr(>F)",i]=anovaD[["Pr(>F)"]][2]
		 # Out-sample
		res_modIout=(t1out-(coefI[1]+insIout%*%coefI[2:length(coefI)]))
		res_modDout=(t1out-(coefD[1]+insDout%*%coefD[2:length(coefD)]))
		anovaMatOut["df_modI",i]=anovaMatIn["df_modI",i]
		anovaMatOut[c("RSS_modD"),i]=var(res_modDout)*(n-1)
		anovaMatOut[c("RSS_modI"),i]=var(res_modIout)*(n-1)
		anovaMatOut["SumOfSqD",i]=(var(res_modDout)-var(res_modIout))*(n-1)
		anovaMatOut["F",i]=((var(res_modDout)-var(res_modIout))*(n-1))/((var(res_modIout)*(n-1))/(n-p-1))
		anovaMatOut["Pr(>F)",i]= 1 - pf(anovaMatOut["F",i],1,anovaMatOut["df_modI",i])
		}

	## Compute in and out of sample f^2
	fsqVec=(rsqVecI-rsqVecD)/(1-rsqVecI) #NOTE!! Using rsqVecLscS (rsq_all) data
	rownames(fsqVec)=c("f.squared","adj.f.squared","out.f.squared","out.adj.f.squared")

	shortList=list()
	shortList[[1]]=rsqVecI
	shortList[[2]]=rsqVecD
	shortList[[3]]=fsqVec
	shortList[[4]]=fsqIdealVec
	shortList[[5]]=anovaMatIn
	shortList[[6]]=anovaMatOut
	shortList[[7]]=resCoefVec
	shortList[[8]]=coefMatD
	names(shortList)=c("rsqVecI","rsqVecD","fsqVecD","fsqIdealVec",
					"anovaMatIn","anovaMatOut","resCoefVec","coefMatD")
	J=jjList[as.character(j)]
	rsqList[[J]]=shortList
	}	
	names(rsqList)=colnames(inS)[jjVec]
	rsqList
	}


##%%%%%%% make cv.rsqList function %%%%%%

makeRsqListCV=function(tar,inS,Imat,ImatOut,jjVec){
	jjList=(1:length(jjVec))
	names(jjList)=as.character(jjVec)
	rsqList=list()

    for(j in jjVec){
	Nmod=ncol(Imat)
	jjj=1:ncol(inS)
	jj=setdiff(jjj,j)
	jjjcoef=1:(ncol(inS)+1)
	jjcoef=setdiff(jjjcoef,j+1)
	coefMatD=matrix(0,length(jjcoef),Nmod)
	rownames(coefMatD)=c("intercept",colnames(inS)[jj])
	colnames(coefMatD)=paste("mod",1:Nmod,sep="")
	rsqVecI=matrix(0,6,Nmod)
	rownames(rsqVecI)=c("r.squared","adj.r.squared","cv.r.squared","cv.adj.r.squared",
					"out.r.squared","out.adj.r.squared")
	rsqVecD=rsqVecI
	fsqIdealVec=matrix(0,1,Nmod)
	rownames(fsqIdealVec)=c("fsqIdeal")

	## residual coefficients are computed on the residuals
	## of the dropped models
	resCoefVec=matrix(0,2,Nmod)
	rownames(resCoefVec)=c("coef","resCoefInD")
	

	## 	n=number of obs; p=number of inputs;
	##	df_modI=(n-p-1); df modD=(n-(p-1)-1)=(n-p)
	##	F-stat  = ((RSS_modD - RSS_modI)/(df_modD=df_modI)) / (RSS_modI/df_modI)
	##          = ((var(res_modD)-var(res_modI))*(n-1)) / ((var(res_modI)*(n-1))/(n-p-1))
	##    Pr(>F)  = 1 - pf(F-stat,df_modD1-df_modS,df_modS)
	
	anovaMatIn=matrix(0,6,Nmod)
	rownames(anovaMatIn)=c("df_modI","RSS_modD","RSS_modI","SumOfSqD","F","Pr(>F)")
	anovaMatOut=anovaMatIn

	for(i in 1:Nmod){
		iiin=Imat[,i]
		t1in=tar[iiin]
		insI=inS[iiin,]
		insD=inS[iiin,jj]
		modD=lm(t1in ~ insD)
		modI=lm(t1in ~insI)
		modDsum=summary(modD)
		modIsum=summary(modI)
		coefD=modD$coefficients
		coefD[!is.finite(coefD)]=0
		coefI=modI$coefficients
		coefI[!is.finite(coefI)]=0
		coefMatD[,i]=coefD
		## rsqVecD[1,i]=modDsum$r.squared
		## rsqVecD[2,i]=modDsum$adj.r.squared
		## rsqVecI[1,i]=modIsum$r.squared
		## rsqVecI[2,i]=modIsum$adj.r.squared
		rsqVecD[1:4,i]=crossValD1(t1in,insD,modD,modDsum)
		rsqVecI[1:4,i]=crossValD1(t1in,insI,modI,modIsum)

		
		## Compute residual coefficient
		resD=modD$residuals
		inD=inS[iiin,j]
		modRes=lm(resD~inD)
		resCoefVec[2,i]=modRes$coefficients[2]
		resCoefVec[1,i]=j2z(coefI[j+1])

		## Compute the out of sample r.squared, r.squared.drop
		iiout=ImatOut[,i] ## Note ImatOut
		t1out=tar[iiout]
		insDout=inS[iiout,jj]
		insIout=inS[iiout,]
		rsqVecD[5,i]=rsqOut(t1out,insDout,coefD)
		rsqVecD[6,i]=adj.rsqOut(t1out,insDout,coefD)
		rsqVecI[5,i]=rsqOut(t1out,insIout,coefI)
		rsqVecI[6,i]=adj.rsqOut(t1out,insIout,coefI)

		
		## Compute fsqIdeal
		v1=inS[iiin,j]
		coef1=coefI[j+1]
		fsqI=(var(v1*coef1)/var(t1in))/(1-rsqVecI[1,i])
		fsqIdealVec[i]=fsqI
	
		## Compute anova SS and F values
		anovaD=anova(modD,modI)
		n=length(iiout);p=sum(j2z(coefI)!=0)
		 # In-sample c("df_modI","RSS_modD","RSS_modI","SumOfSqD","F","Pr(>F)")
		anovaMatIn["df_modI",i]=anovaD[["Res.Df"]][2]
		anovaMatIn[c("RSS_modD","RSS_modI"),i]=anovaD[["RSS"]]
		anovaMatIn["SumOfSqD",i]=anovaD[["Sum of Sq"]][2]
		anovaMatIn["F",i]=anovaD[["F"]][2]
		anovaMatIn["Pr(>F)",i]=anovaD[["Pr(>F)"]][2]
		 # Out-sample
		res_modIout=(t1out-(coefI[1]+insIout%*%coefI[2:length(coefI)]))
		res_modDout=(t1out-(coefD[1]+insDout%*%coefD[2:length(coefD)]))
		anovaMatOut["df_modI",i]=anovaMatIn["df_modI",i]
		anovaMatOut[c("RSS_modD"),i]=var(res_modDout)*(n-1)
		anovaMatOut[c("RSS_modI"),i]=var(res_modIout)*(n-1)
		anovaMatOut["SumOfSqD",i]=(var(res_modDout)-var(res_modIout))*(n-1)
		anovaMatOut["F",i]=((var(res_modDout)-var(res_modIout))*(n-1))/((var(res_modIout)*(n-1))/(n-p-1))
		anovaMatOut["Pr(>F)",i]= 1 - pf(anovaMatOut["F",i],1,anovaMatOut["df_modI",i])
		}

	## Compute in and out of sample f^2
	fsqVec=(rsqVecI-rsqVecD)/(1-rsqVecI) #NOTE!! Using rsqVecLscS (rsq_all) data
	rownames(fsqVec)=c("f.squared","adj.f.squared","cv.f.squared","cv.adj.f.squared","out.f.squared","out.adj.f.squared")

	shortList=list()
	shortList[[1]]=rsqVecI
	shortList[[2]]=rsqVecD
	shortList[[3]]=fsqVec
	shortList[[4]]=fsqIdealVec
	shortList[[5]]=anovaMatIn
	shortList[[6]]=anovaMatOut
	shortList[[7]]=resCoefVec
	shortList[[8]]=coefMatD
	names(shortList)=c("rsqVecI","rsqVecD","fsqVecD","fsqIdealVec",
					"anovaMatIn","anovaMatOut","resCoefVec","coefMatD")
	J=jjList[as.character(j)]
	rsqList[[J]]=shortList
	}	
	names(rsqList)=colnames(inS)[jjVec]
	rsqList
	}


##%%%%%%%%%%%% rsq data computations compiled into matrices %%%%%


##%%%%%%%%%%%%%% Make rsq Matrices and Array's %%%%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## Function to insure that all inputs are good
nonRed=function(taR,inputS){
	mod=lm(taR~inputS)
	coeF=mod$coef
	goodIns=inputS[,is.finite(coeF[2:length(coeF)])]
	goodIns
	}

## Function to compute r^2, f^2, for in, out, outL, adj,
## cv, cv10, cv5, cv2 for all inputs with p-value<0.1

makeRsqMat0.1=function(tar,inS,IIin,IIout,IIoutL){
	goodIns=colnames(nonRed(tar[IIin],inS[IIin,]))
	inS=inS[,goodIns]
	modI=lm(tar[IIin]~inS[IIin,])
	sumModI=summary(modI)
	n=length(tar[IIin])
	p=nrow(sumModI$coef)-1
	df_modI=n-p-1
	df_modD=n-p	
	## 	n=number of obs; p=number of inputs;
	##	df_modI=(n-p-1); df modD=(n-(p-1)-1)=(n-p)
	##	F-stat  = ((RSS_modD - RSS_modI)/(df_modD-df_modI)) / (RSS_modI/df_modI)
	##          = ((var(res_modD)-var(res_modI))*(n-1)) / ((var(res_modI)*(n-1))/(n-p-1))
	##          =df_modI * f.squared
	##    Pr(>F)  = 1 - pf(F-stat,df_modD-df_modI,df_modI)
	##            = 1 - pf(F-stat,1,df_modI)
	inVec=c(n,p,df_modI,sumModI$r.squared,sumModI$adj.r.squared)
	inCoef=modI$coef
	inCoef[!is.finite(inCoef)]=0
	
	## Compute out of sample r^2
	out.r.squaredI=rsqOut(tar[IIout],inS[IIout,],inCoef)
	outL.r.squaredI=rsqOut(tar[IIoutL],inS[IIoutL,],inCoef)
	inOutVec=c(inVec,out.r.squaredI,outL.r.squaredI)
	names(inOutVec)=c("n-obs","n-inputs","df_modI","r.squaredI",
					"adj.r.squaredI","out.r.squaredI","outL.r.squaredI")
	
	## Find inputs with p-value <= 0.1
	sumCoef=(sumModI$coef)[2:nrow(sumModI$coef),]
	colNamesSum=rownames(sumModI$coef)[2:nrow(sumModI$coef)]
	colNamesSum=substring(colNamesSum,nchar("inS[IIin, ]")+1,nchar(colNamesSum))
	rownames(sumCoef)=colNamesSum
	jjVec0.1=(1:nrow(sumCoef))[sumCoef[,4]<=0.1]
	names(jjVec0.1)=colNamesSum[jjVec0.1]
	
	## Compute in/out r^2 for dropped variables
	colNamesD=c("df_modD","r.squaredD","adj.r.squaredD","out.r.squaredD",
				"outL.r.squaredD",colnames(sumCoef))
	rowNamesD=names(jjVec0.1)
	inOutMatD=matrix(0,length(jjVec0.1),length(colNamesD))
	rownames(inOutMatD)=rowNamesD
	colnames(inOutMatD)=colNamesD
	
	for(s in 1:length(jjVec0.1)){
		j=jjVec0.1[s]
		tVec=sumCoef[j,]
		jjD=setdiff(1:ncol(inS),j)
		modD=lm(tar[IIin]~inS[IIin,jjD])
		sumModD=summary(modD)
		## Collect basic model stats
		inVecD=c(df_modD,sumModD$r.squared,sumModD$adj.r.squared)
		inCoefD=modD$coef
		inCoefD[!is.finite(inCoefD)]=0
	
		## Compute out of sample r^2
		out.r.squaredD=rsqOut(tar[IIout],inS[IIout,jjD],inCoefD)
		outL.r.squaredD=rsqOut(tar[IIoutL],inS[IIoutL,jjD],inCoefD)
		inOutVecD=c(inVecD,out.r.squaredD,outL.r.squaredD)
		inOutMatD[s,1:length(inOutVecD)]=inOutVecD
		inOutMatD[s,(length(inOutVecD)+1):ncol(inOutMatD)]=tVec
		}
	
	## Compute cross validation values
	cv.rsq=t(crossValD(tar[IIin],inS[IIin,],jjVec0.1))
	cv10.rsq=t(crossValM(tar[IIin],inS[IIin,],jjVec0.1,10))
	cv5.rsq=t(crossValM(tar[IIin],inS[IIin,],jjVec0.1,5))
	cv2.rsq=t(crossValM(tar[IIin],inS[IIin,],jjVec0.1,2))
	colnames(cv10.rsq)=c("r.squaredI","r.squaredD","cv10.r.squaredI",
				    "cv10.r.squaredD","f.squared","cv10.f.squared")
	colnames(cv5.rsq)=c("r.squaredI","r.squaredD","cv5.r.squaredI",
				    "cv5.r.squaredD","f.squared","cv5.f.squared")
	colnames(cv2.rsq)=c("r.squaredI","r.squaredD","cv2.r.squaredI",
				    "cv2.r.squaredD","f.squared","cv2.f.squared")

	##%% inOutVec to matrix
	inOutMat=matrix(1,length(jjVec0.1),1)%*%inOutVec
	rownames(inOutMat)=names(jjVec0.1)
	colnames(inOutMat)=names(inOutVec)
	
	## collect values into common matrices
	## rsqI
	rsqMatI=cbind(inOutMat,cv.rsq[,"cv.r.squaredI"],cv10.rsq[,"cv10.r.squaredI"],
		cv5.rsq[,"cv5.r.squaredI"],cv2.rsq[,"cv2.r.squaredI"])
	colnames(rsqMatI)[(ncol(rsqMatI)-3):ncol(rsqMatI)]=c("cv.r.squaredI",
					"cv10.r.squaredI","cv5.r.squaredI","cv2.r.squaredI")
	rsqMatI=cbind(matrix(length(jjVec0.1),nrow(rsqMatI),1),rsqMatI)
	colnames(rsqMatI)[1]="numb. inputs w. pval<=0.1"

	## rsqD
	rsqMatD=cbind(inOutMatD[,c(1,6:9,2:5)],cv.rsq[,"cv.r.squaredD"],
		cv10.rsq[,"cv10.r.squaredD"],cv5.rsq[,"cv5.r.squaredD"],cv2.rsq[,"cv2.r.squaredD"])
	colnames(rsqMatD)[(ncol(rsqMatD)-3):ncol(rsqMatD)]=c("cv.r.squaredD",
					"cv10.r.squaredD","cv5.r.squaredD","cv2.r.squaredD")
	## f.squared
	rsqIcols=c("r.squaredI","adj.r.squaredI","out.r.squaredI","outL.r.squaredI","cv.r.squaredI",  
			"cv10.r.squaredI","cv5.r.squaredI","cv2.r.squaredI")
	rsqDcols=c("r.squaredD","adj.r.squaredD","out.r.squaredD","outL.r.squaredD","cv.r.squaredD",  
			"cv10.r.squaredD","cv5.r.squaredD","cv2.r.squaredD")
	fsqMat=(rsqMatI[,rsqIcols]-rsqMatD[,rsqDcols])/(1-rsqMatI[,rsqIcols])
	colnames(fsqMat)=c("f.squared","adj.f.squared","out.f.squared","outL.f.squared","cv.f.squared",  
			"cv10.f.squared","cv5.f.squared","cv2.f.squared")
	
	## Fstat
	dfMat=rsqMatI[,"df_modI"]%*%matrix(1,1,ncol(fsqMat))
	fstatMat=dfMat*fsqMat
	colnames(fstatMat)=c("fstat","adj.fstat","out.fstat","outL.fstat","cv.fstat",  
			"cv10.fstat","cv5.fstat","cv2.fstat")

	## Fstat p-values
	fpvalMat=1-pf(fstatMat,1,rsqMatI[1,"df_modI"])
	colnames(fpvalMat)=c("fstat p-value","adj.fstat p-value","out.fstat p-value","outL.fstat p-value",
				"cv.fstat p-value","cv10.fstat p-value","cv5.fstat p-value","cv2.fstat p-value")


	## FtoT
	tstatFmat=FtoT(fstatMat,1,rsqMatI[1,"df_modI"])
	colnames(tstatFmat)=c("tstatF","adj.tstatF","out.tstatF","outL.tstatF","cv.tstatF",  
			"cv10.tstatF","cv5.tstatF","cv2.tstatF")

	## fstat - lower bounds
	lowb05=qf(0.05,1,rsqMatI[1,"df_modI"],fstatMat[,"fstat"])
	lowb05[lowb05>fstatMat[,"fstat"]]=0 # in case lowb>fstat stat due to numerical issues
	lowb01=qf(0.01,1,rsqMatI[1,"df_modI"],fstatMat[,"fstat"])
	lowb01[lowb01>fstatMat[,"fstat"]]=0 # in case lowb>fstat due to numerical issues
	lowbMat=cbind(lowb05,lowb01)
	colnames(lowbMat)=c("in.fstat 95% lower bound","in.fstat 99% lower bound")
	rownames(lowbMat)=rownames(rsqMatI)

	## prob out, outL less than in.fstat lower bound, less than zero
	confIntMat=cbind(fstatMat[,"out.fstat"]<lowbMat[,1],fstatMat[,"out.fstat"]<lowbMat[,2],
				 fstatMat[,"out.fstat"]<=0,fstatMat[,"outL.fstat"]<lowbMat[,1],
					fstatMat[,"outL.fstat"]<lowbMat[,2],fstatMat[,"outL.fstat"]<=0)
	colnames(confIntMat)=c("pr out.fstat<95% lowb","pr out.fstat<99% lowb","pr out.fstat<=0",
					"pr outL.fstat<95% lowb","pr outL.fstat<99% lowb","pr outL.fstat<=0")
	rownames(confIntMat)=rownames(rsqMatI)
	
	##%%%%%%%%%% Put all this together %%%%
	outMatA=cbind(rsqMatI,rsqMatD,fsqMat,fstatMat,fpvalMat,tstatFmat,lowbMat,confIntMat)
	outMatA
	}




## Function to do these comutations over a large set of models

makeRsqMat0.1L=function(tar,inS,Imat,ImatOut,ImatOutL){
	countS1=paste("mod00",1:9,sep="")
	countS2=paste("mod0",10:99,sep="")
	countS3=paste("mod",100:999,sep="")
	countSA=c(countS1,countS2,countS3)
	nMod=ncol(Imat)
	modNames=countSA[1:nMod]
	
	outMatS=makeRsqMat0.1(tar,inS,Imat[,1],ImatOut[,1],ImatOutL[,1])
	rownames(outMatS)=paste(modNames[1],rownames(outMatS))
	outMat=outMatS
	for(i in 2:nMod){
		## Test for (npval<0.1)<2
		sumModCoef=summary(lm(tar[Imat[,i]]~inS[Imat[,i],]))$coef
		sumModCoef=sumModCoef[2:nrow(sumModCoef),]
		if(sum(sumModCoef[,4]<=0.1)<2){i=i+1}
		flush.console()
		print(i)
		outMatS=makeRsqMat0.1(tar,inS,Imat[,i],ImatOut[,i],ImatOutL[,i])
		rownames(outMatS)=paste(modNames[i],rownames(outMatS))
		outMat=rbind(outMat,outMatS)
		}
	outMat
	}


## Function to compute r^2, f^2, for in, out, outL, adj, cv
## for all inputs

makeRsqMatNcvA=function(tar,inS,IIin,IIout,IIoutL){
	goodIns=colnames(nonRed(tar[IIin],inS[IIin,]))
	inS=inS[,goodIns]
	modI=lm(tar[IIin]~inS[IIin,])
	sumModI=summary(modI)
	n=length(tar[IIin])
	p=nrow(sumModI$coef)-1
	df_modI=n-p-1
	df_modD=n-p	
	## 	n=number of obs; p=number of inputs;
	##	df_modI=(n-p-1); df modD=(n-(p-1)-1)=(n-p)
	##	F-stat  = ((RSS_modD - RSS_modI)/(df_modD-df_modI)) / (RSS_modI/df_modI)
	##          = ((var(res_modD)-var(res_modI))*(n-1)) / ((var(res_modI)*(n-1))/(n-p-1))
	##          =df_modI * f.squared
	##    Pr(>F)  = 1 - pf(F-stat,df_modD-df_modI,df_modI)
	##            = 1 - pf(F-stat,1,df_modI)
	inVec=c(n,p,df_modI,sumModI$r.squared,sumModI$adj.r.squared)
	inCoef=modI$coef
	inCoef[!is.finite(inCoef)]=0
	
	## Compute out of sample r^2
	out.r.squaredI=rsqOut(tar[IIout],inS[IIout,],inCoef)
	outL.r.squaredI=rsqOut(tar[IIoutL],inS[IIoutL,],inCoef)
	inOutVec=c(inVec,out.r.squaredI,outL.r.squaredI)
	names(inOutVec)=c("n-obs","n-inputs","df_modI","r.squaredI",
					"adj.r.squaredI","out.r.squaredI","outL.r.squaredI")
	
	## Find inputs with p-value <= 0.1
	sumCoef=(sumModI$coef)[2:nrow(sumModI$coef),]
	colNamesSum=rownames(sumModI$coef)[2:nrow(sumModI$coef)]
	colNamesSum=substring(colNamesSum,nchar("inS[IIin, ]")+1,nchar(colNamesSum))
	rownames(sumCoef)=colNamesSum
	nPval=sum(sumCoef[,4]<=0.1)

	## Find inputs for goodIns
	jjVec0.1=1:ncol(inS)
	names(jjVec0.1)=colnames(inS)
	
	## Compute in/out r^2 for dropped variables
	colNamesD=c("df_modD","r.squaredD","adj.r.squaredD","out.r.squaredD",
				"outL.r.squaredD",colnames(sumCoef))
	rowNamesD=names(jjVec0.1)
	inOutMatD=matrix(0,length(jjVec0.1),length(colNamesD))
	rownames(inOutMatD)=rowNamesD
	colnames(inOutMatD)=colNamesD
	
	for(s in 1:length(jjVec0.1)){
		j=jjVec0.1[s]
		tVec=sumCoef[j,]
		jjD=setdiff(1:ncol(inS),j)
		modD=lm(tar[IIin]~inS[IIin,jjD])
		sumModD=summary(modD)
		## Collect basic model stats
		inVecD=c(df_modD,sumModD$r.squared,sumModD$adj.r.squared)
		inCoefD=modD$coef
		inCoefD[!is.finite(inCoefD)]=0
	
		## Compute out of sample r^2
		out.r.squaredD=rsqOut(tar[IIout],inS[IIout,jjD],inCoefD)
		outL.r.squaredD=rsqOut(tar[IIoutL],inS[IIoutL,jjD],inCoefD)
		inOutVecD=c(inVecD,out.r.squaredD,outL.r.squaredD)
		inOutMatD[s,1:length(inOutVecD)]=inOutVecD
		inOutMatD[s,(length(inOutVecD)+1):ncol(inOutMatD)]=tVec
		}
	
	## Compute cross validation values
	cv.rsq=t(crossValD(tar[IIin],inS[IIin,],jjVec0.1))

	##%% inOutVec to matrix
	inOutMat=matrix(1,length(jjVec0.1),1)%*%inOutVec
	rownames(inOutMat)=names(jjVec0.1)
	colnames(inOutMat)=names(inOutVec)
	
	## collect values into common matrices
	## rsqI
	rsqMatI=cbind(inOutMat,cv.rsq[,"cv.r.squaredI"])
	colnames(rsqMatI)[ncol(rsqMatI)]=c("cv.r.squaredI")
	rsqMatI=cbind(matrix(length(jjVec0.1),nrow(rsqMatI),1),rsqMatI)
	colnames(rsqMatI)[1]="numb. inputs w. pval<=0.1"
	rsqMatI[,1]=nPval

	## rsqD
	rsqMatD=cbind(inOutMatD[,c(1,6:9,2:5)],cv.rsq[,"cv.r.squaredD"])
	colnames(rsqMatD)[ncol(rsqMatD)]=c("cv.r.squaredD")
	## f.squared
	rsqIcols=c("r.squaredI","adj.r.squaredI","out.r.squaredI","outL.r.squaredI","cv.r.squaredI")
	rsqDcols=c("r.squaredD","adj.r.squaredD","out.r.squaredD","outL.r.squaredD","cv.r.squaredD")
	fsqMat=(rsqMatI[,rsqIcols]-rsqMatD[,rsqDcols])/(1-rsqMatI[,rsqIcols])
	colnames(fsqMat)=c("f.squared","adj.f.squared","out.f.squared","outL.f.squared","cv.f.squared")
	
	## Fstat
	dfMat=rsqMatI[,"df_modI"]%*%matrix(1,1,ncol(fsqMat))
	fstatMat=dfMat*fsqMat
	colnames(fstatMat)=c("fstat","adj.fstat","out.fstat","outL.fstat","cv.fstat")

	## Fstat p-values
	fpvalMat=1-pf(fstatMat,1,rsqMatI[1,"df_modI"])
	colnames(fpvalMat)=c("fstat p-value","adj.fstat p-value","out.fstat p-value","outL.fstat p-value",
				"cv.fstat p-value")


	## FtoT
	tstatFmat=FtoT(fstatMat,1,rsqMatI[1,"df_modI"])
	colnames(tstatFmat)=c("tstatF","adj.tstatF","out.tstatF","outL.tstatF","cv.tstatF")

	## fstat - lower bounds
	lowb05=qf(0.05,1,rsqMatI[1,"df_modI"],fstatMat[,"fstat"])
	lowb05[lowb05>fstatMat[,"fstat"]]=0 # in case lowb>fstat stat due to numerical issues
	lowb01=qf(0.01,1,rsqMatI[1,"df_modI"],fstatMat[,"fstat"])
	lowb01[lowb01>fstatMat[,"fstat"]]=0 # in case lowb>fstat due to numerical issues
	lowbMat=cbind(lowb05,lowb01)
	colnames(lowbMat)=c("in.fstat 95% lower bound","in.fstat 99% lower bound")
	rownames(lowbMat)=rownames(rsqMatI)

	## prob out, outL less than in.fstat lower bound, less than zero
	confIntMat=cbind(fstatMat[,"out.fstat"]<lowbMat[,1],fstatMat[,"out.fstat"]<lowbMat[,2],
				 fstatMat[,"out.fstat"]<=0,fstatMat[,"outL.fstat"]<lowbMat[,1],
					fstatMat[,"outL.fstat"]<lowbMat[,2],fstatMat[,"outL.fstat"]<=0)
	colnames(confIntMat)=c("pr out.fstat<95% lowb","pr out.fstat<99% lowb","pr out.fstat<=0",
					"pr outL.fstat<95% lowb","pr outL.fstat<99% lowb","pr outL.fstat<=0")
	rownames(confIntMat)=rownames(rsqMatI)
	
	##%%%%%%%%%% Put all this together %%%%
	outMatA=cbind(rsqMatI,rsqMatD,fsqMat,fstatMat,fpvalMat,tstatFmat,lowbMat,confIntMat)
	outMatA
	}


## Function to do these comutations over a large set of models

makeRsqMatNcvAL=function(tar,inS,Imat,ImatOut,ImatOutL){
	countS1=paste("mod00",1:9,sep="")
	countS2=paste("mod0",10:99,sep="")
	countS3=paste("mod",100:999,sep="")
	countSA=c(countS1,countS2,countS3)
	nMod=ncol(Imat)
	modNames=countSA[1:nMod]
	
	outMatS=makeRsqMatNcvA(tar,inS,Imat[,1],ImatOut[,1],ImatOutL[,1])
	rownames(outMatS)=paste(modNames[1],rownames(outMatS))
	outMat=outMatS
	for(i in 2:nMod){
		## Test for (npval<0.1)<2
		sumModCoef=summary(lm(tar[Imat[,i]]~inS[Imat[,i],]))$coef
		sumModCoef=sumModCoef[2:nrow(sumModCoef),]
		if(sum(sumModCoef[,4]<=0.1)<2){i=i+1}
		flush.console()
		print(i)
		outMatS=makeRsqMatNcvA(tar,inS,Imat[,i],ImatOut[,i],ImatOutL[,i])
		rownames(outMatS)=paste(modNames[i],rownames(outMatS))
		outMat=rbind(outMat,outMatS)
		}
	outMat
	}





##%%%%%%%%%% Fuction to compute fstat mats from rsqMat's %%%%%%%%%%%%

makeFstatMat=function(mbOut,DF,prefix="in."){
	foo=mbOut
	junk1=round(cbind(foo[[2]],foo[[1]],t(foo[[3]])),3)
	junk2=round(1-pf(cbind(foo[[1]],t(foo[[3]])),1,DF),4)
	junk3=round(FtoT(cbind(foo[[1]],t(foo[[3]])),1,DF),2)
	junkA=cbind(junk1,junk2,junk3)
	colnames(junkA)=c("std outL.fstat","outL.fstat",paste(prefix,"fstat",sep=""),"outL.pval",
			paste(prefix,"pval",sep=""),"outL.tstatF",paste(prefix,"tstatF",sep=""))
	junkA
	}


##%%%%%%% function to collect stats on individual inputs per sim %%%%%%%%


inpStats=function(rsqLS,inpN,funC=Cmedian){
	inpNU=unique(inpN)
	outMat=matrix(0,length(inpNU),ncol(rsqLS)+1)
	for(i in 1:length(inpNU)){
		inp1=inpNU[i]
		jj=(inpN==inp1)
		outMat[i,1]=sum(jj)
		rsqLSS=rsqLS[jj,]
		outMat[i,2:(ncol(outMat)-6)]=funC(rsqLSS[,1:(ncol(rsqLSS)-6)])
		outMat[i,(ncol(outMat)-5):ncol(outMat)]=Cmean(rsqLSS[,(ncol(rsqLSS)-5):ncol(rsqLSS)])
		}
	rownames(outMat)=inpNU
	colnames(outMat)=c("n-pval<=0.1",colnames(rsqLS))
	outMat
	}

inpStatsFunc=function(rsqLS,inpN,funC=Cmedian){
	inpNU=unique(inpN)
	outMat=matrix(0,length(inpNU),ncol(rsqLS)+1)
	for(i in 1:length(inpNU)){
		inp1=inpNU[i]
		jj=(inpN==inp1)&(rsqLS[,"Pr(>|t|)"]<=0.1)
		outMat[i,1]=sum(jj)
		rsqLSS=rsqLS[jj,]
		outMat[i,2:(ncol(outMat)-6)]=funC(rsqLSS[,1:(ncol(rsqLSS)-6)])
		outMat[i,(ncol(outMat)-5):ncol(outMat)]=Cmean(rsqLSS[,(ncol(rsqLSS)-5):ncol(rsqLSS)])
		}
	rownames(outMat)=inpNU
	colnames(outMat)=c("n-pval<=0.1",colnames(rsqLS))
	outMat
	}


##%%%%%%%%%%% Function for Analysis of rsqList's %%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

makeAnalysisMat=function(rsqList,jjVec){
	outMat=matrix(0,15,length(rsqList))
	for(i in 1:length(rsqList)){
		foo0=Rmedian(rsqList[[i]][[1]])
		ii=rsqList[[i]][[5]][6,]<0.05
		foo1=as.matrix(median(rsqList[[i]][[7]][2,]/rsqList[[i]][[7]][1,]))
		rownames(foo1)=c("median-(resCoefInD/in.coefInD)")
		foo1.1=as.matrix(median(rsqList[[i]][[7]][1,]))
		rownames(foo1.1)=c("median-in.coefInD")
		foo2=as.matrix(sum(ii)/length(ii))
		rownames(foo2)=c("count[ii]/count[mods], ii=mods-Pr(>in.F)<0.05")
		foo3=Rmedian(rsqList[[i]][[3]][,ii])
		rownames(foo3)=paste("median-",rownames(foo3),"-on Pr(>in.F)<0.05",sep="")
		foo4=as.matrix(sum(rsqList[[i]][[6]][5,ii]<0)/sum(ii))
		rownames(foo4)=c("percent-out.Fstat<0 on Pr(>in.F)<0.05")
		foo5=Rsum(rsqList[[i]][[3]][2:4,ii]<0)/sum(ii)
		rownames(foo5)=paste("percent-",rownames(foo5),"<0-on Pr(>in.F)<0.05",sep="")
		fooOut=rbind(foo0,foo1,foo1.1,foo2,foo3,foo4,foo5)
		outMat[,i]=fooOut
		}
	rownames(outMat)=rownames(fooOut)
	colnames(outMat)=names(jjVec)
	outMat
	}

## A similar analysis function, but with a focus on out.f^2<0 
## for both Pr(>in.F)<0.05 and Pr(>in.F)<0.01




makeAnalysisMatI=function(rsqList,jjVec){
	outMat=matrix(0,9,length(rsqList))
	for(i in 1:length(rsqList)){
		# foo0=Rmedian(rsqList[[i]][[1]])
		ii=rsqList[[i]][[5]][6,]<0.05
		jj=rsqList[[i]][[5]][6,]<0.01
		# foo1=as.matrix(median(rsqList[[i]][[7]][2,]/rsqList[[i]][[7]][1,]))
		# rownames(foo1)=c("median-(resCoefInD/in.coefInD)")
		foo1=as.matrix(median(rsqList[[i]][[5]][6,]))
		rownames(foo1)=c("median-Pr(>in.F)")
		# foo1.1=as.matrix(median(rsqList[[i]][[7]][1,]))
		# rownames(foo1.1)=c("median-in.coefInD")
		foo2=as.matrix(sum(ii)/length(ii))
		rownames(foo2)=c("count[ii]/count[mods], ii=Pr(>in.F)<0.05")
		foo3=Rmedian(rsqList[[i]][[3]][c(1,3),ii])
		rownames(foo3)=paste("median-",rownames(foo3),"-on Pr(>in.F)<0.05",sep="")
		# foo4=as.matrix(sum(rsqList[[i]][[6]][5,ii]<0)/sum(ii))
		# rownames(foo4)=c("percent-out.Fstat<0 on Pr(>in.F)<0.05")
		foo5=Rsum(rsqList[[i]][[3]][3:4,ii]<0)/sum(ii)
		rownames(foo5)=paste("percent-",rownames(foo5),"<0-on Pr(>in.F)<0.05",sep="")
		
		## Start jj section
		foo6=as.matrix(sum(jj)/length(jj))
		rownames(foo6)=c("count[jj]/count[mods], jj=Pr(>in.F)<0.01")
		foo7=Rmedian(rsqList[[i]][[3]][c(1,3),jj])
		rownames(foo7)=paste("median-",rownames(foo7),"-on Pr(>in.F)<0.01",sep="")
		foo8=Rsum(rsqList[[i]][[3]][3:4,jj]<0)/sum(jj)
		rownames(foo8)=paste("percent-",rownames(foo8),"<0-on Pr(>in.F)<0.01",sep="")
		fooOut=rbind(foo1,foo2,foo3,foo5,foo6,foo7,foo8)
		outMat[,i]=fooOut[c(1:5,7:10),]
		}
	rownames(outMat)=rownames(fooOut)[c(1:5,7:10)]
	colnames(outMat)=names(jjVec)
	outMat
	}


##%%%%%%%%%%%%%%%% Make Confidence interval matrix %%%%%%%%%%%%%
## Function to compute two sided confidence intervals

makeConfIntMat=function(rsqList){
	confIntMat=matrix(0,5,length(rsqList))
      rownames(confIntMat)=c("median in.Pr(F>0)","pr out.F outside 90% confInt",
					"pr out.F outside 95% confInt","pr out.F outside 99% confInt",
					"pr out.F<0")

	for(i in 1:length(rsqList)){
		medFstatProb=median(rsqList[[i]][[5]][6,])
		in.fstat=rsqList[[i]][[5]][5,]
		out.fstat=rsqList[[i]][[6]][5,]
		df=rsqList[[i]][[5]]["df_modI",1]
		lowb1=qf(0.05,1,df,in.fstat)
		lowb05=qf(0.025,1,df,in.fstat)
		lowb01=qf(0.005,1,df,in.fstat)
		upb1=qf(0.95,1,df,in.fstat)
		upb05=qf(0.975,1,df,in.fstat)
		upb01=qf(0.995,1,df,in.fstat)

		confIntMat[1,i]=medFstatProb
		confIntMat[2,i]=sum((out.fstat<lowb1)|(out.fstat>upb1))/length(out.fstat)
		confIntMat[3,i]=sum((out.fstat<lowb05)|(out.fstat>upb05))/length(out.fstat)
		confIntMat[4,i]=sum((out.fstat<lowb01)|(out.fstat>upb01))/length(out.fstat)
		confIntMat[5,i]=sum(out.fstat < 0)/length(out.fstat)
		}
	colnames(confIntMat)=names(rsqList)
	confIntMat
	}


## Function to coverage for lower bound confidence intervals

makeConfIntMatL=function(rsqList){
	confIntMat=matrix(0,5,length(rsqList))
      rownames(confIntMat)=c("median in.Pr(F>0)","pr out.F below 90% lowb",
					"pr out.F below 95% lowb","pr out.F below 99% lowb",
					"pr out.F<0")

	for(i in 1:length(rsqList)){
		medFstatProb=median(rsqList[[i]][[5]][6,])
		in.fstat=rsqList[[i]][[5]][5,]
		out.fstat=rsqList[[i]][[6]][5,]
		df=rsqList[[i]][[5]]["df_modI",1]
		lowb1=qf(0.1,1,df,in.fstat)
		lowb05=qf(0.05,1,df,in.fstat)
		lowb01=qf(0.01,1,df,in.fstat)

		confIntMat[1,i]=medFstatProb
		confIntMat[2,i]=sum(out.fstat < lowb1)/length(out.fstat)
		confIntMat[3,i]=sum(out.fstat < lowb05)/length(out.fstat)
		confIntMat[4,i]=sum(out.fstat < lowb01)/length(out.fstat)
		confIntMat[5,i]=sum(out.fstat < 0)/length(out.fstat)
		}
	colnames(confIntMat)=names(rsqList)
	confIntMat
	}


##%%%% Function to compute coverage prob. just where Pr(>F)<=0.05

makeConfIntMatL05=function(rsqList){
	confIntMat=matrix(0,6,length(rsqList))
      rownames(confIntMat)=c("median in.Pr(F)","pr in.Pr(F)<0.05, =ii05",
					"pr out.F below 95% lowb ii05","pr out.F below 99% lowb ii05",
					"pr out.F<=0 ii05","median(out.f^2)/median(in.f^2) ii05")

	for(i in 1:length(rsqList)){
		medFstatProb=median(rsqList[[i]][[5]][6,])
		ii05=rsqList[[i]][[5]][6,]<=0.05
		in.fstat=rsqList[[i]][[5]][5,]
		out.fstat=rsqList[[i]][[6]][5,]
		df=rsqList[[i]][[5]]["df_modI",1]
		lowb1=qf(0.1,1,df,in.fstat)
		lowb05=qf(0.05,1,df,in.fstat)
		lowb01=qf(0.01,1,df,in.fstat)
		lowb00=0
		confIntMat[1,i]=medFstatProb
		confIntMat[2,i]=sum(ii05)/length(ii05)
		confIntMat[3,i]=sum(out.fstat[ii05] < lowb05[ii05])/sum(ii05)
		confIntMat[4,i]=sum(out.fstat[ii05] < lowb01[ii05])/sum(ii05)
		confIntMat[5,i]=sum(out.fstat[ii05] <= 0)/sum(ii05)
		confIntMat[6,i]=median(rsqList[[i]][[3]][3,ii05])/median(rsqList[[i]][[3]][1,ii05])
		}
	colnames(confIntMat)=names(rsqList)
	confIntMat
	}



##%%%%%%%%% Function to compute conf. intervals with cv-data %%%%%

makeConfIntMatLCV=function(rsqList,cv.f){
	confIntMat=matrix(0,6,length(rsqList))
      rownames(confIntMat)=c("median in.Pr(F>0)","pr in.Pr(F)<0.05, =ii05",
					"pr out.F below 95% lowb",
					"pr cv.F below 95% lowb","pr out.F<0",
					"pr cv.F<0")

	for(i in 1:length(rsqList)){
		medFstatProb=median(rsqList[[i]][[5]][6,])
		ii05=rsqList[[i]][[5]][6,]<=0.05
		in.fstat=rsqList[[i]][[5]][5,]
		out.fstat=rsqList[[i]][[6]][5,]
		df=rsqList[[i]][[5]]["df_modI",1]
		cv.fstat=rsqList[[i]][[3]][cv.f,]*df
		lowb05=qf(0.05,1,df,in.fstat)

		confIntMat[1,i]=medFstatProb
		confIntMat[2,i]=sum(ii05)/length(out.fstat)
		confIntMat[3,i]=sum(out.fstat < lowb05)/length(out.fstat)
		confIntMat[4,i]=sum(cv.fstat < lowb05)/length(out.fstat)
		confIntMat[5,i]=sum(out.fstat < 0)/length(out.fstat)
		confIntMat[6,i]=sum(cv.fstat < 0)/length(out.fstat)
		}
	colnames(confIntMat)=names(rsqList)
	confIntMat
	}



makeConfIntMatLCV05=function(rsqList,cv.f){
	confIntMat=matrix(0,6,length(rsqList))
      rownames(confIntMat)=c("median in.Pr(F)","pr in.Pr(F)<0.05, =ii05",
					"pr out.Fstat below 95% lowb ii05","pr cv.Fstat below 95% lowb ii05",
					"pr out.Fstat<0 ii05","pr cv.Fstat<0 ii05")

	for(i in 1:length(rsqList)){
		medFstatProb=median(rsqList[[i]][[5]][6,])
		ii05=rsqList[[i]][[5]][6,]<=0.05
		in.fstat=rsqList[[i]][[5]][5,]
		out.fstat=rsqList[[i]][[6]][5,]
		df=rsqList[[i]][[5]]["df_modI",1]
		cv.fstat=rsqList[[i]][[3]][cv.f,]*df
		lowb05=qf(0.05,1,df,in.fstat)

		confIntMat[1,i]=medFstatProb
		confIntMat[2,i]=sum(ii05)/length(ii05)
		confIntMat[3,i]=sum(out.fstat[ii05] < lowb05[ii05])/sum(ii05)
		confIntMat[4,i]=sum(cv.fstat[ii05] < lowb05[ii05])/sum(ii05)
		confIntMat[5,i]=sum(out.fstat[ii05] < 0)/sum(ii05)
		confIntMat[6,i]=sum(cv.fstat[ii05] < 0)/sum(ii05)
		}
	colnames(confIntMat)=names(rsqList)
	confIntMat
	}

##%%%%%%% Function to compute coverage probability for all cv,outL

makeConfIntMatLCVA=function(cv.rsqList){
	nCvOut=nrow(cv.rsqList[[1]][[3]])
	CvOutNames=rownames(cv.rsqList[[1]][[3]])
	rsqNames=rownames(cv.rsqList[[1]][[1]])
 	confIntMat=matrix(0,3+(7*nCvOut),length(cv.rsqList))
      	 rownames(confIntMat)=c("median in.Pr(F>0)","pr in.Pr(F)<0.05 =ii05",
			"degrees of freedom",
			paste("median",rsqNames),
			paste("median",CvOutNames),
			paste("median tstat of",CvOutNames),
 			paste("pr",CvOutNames,"below 95% lowb"),
 			paste("pr",CvOutNames,"<0"),
 			paste("pr",CvOutNames,"below 95% lowb ii05"),
 			paste("pr",CvOutNames,"<0 ii05")
			)
 
 	for(i in 1:length(cv.rsqList)){
		medFstatProb=median(cv.rsqList[[i]][[5]][6,])
 		ii05=cv.rsqList[[i]][[5]][6,]<=0.05
 		in.fstat=cv.rsqList[[i]][[5]][5,]
 		df=cv.rsqList[[i]][[5]]["df_modI",1]
 		other.fstat=cv.rsqList[[i]][[3]]*df
		lowb05=qf(0.05,1,df,in.fstat)
		lowb05[in.fstat<=lowb05]=0.0001
 		lowb05M=matrix(1,nCvOut,1)%*%lowb05
 
 		confIntMat[1,i]=medFstatProb
 		confIntMat[2,i]=sum(ii05)/length(ii05)
		confIntMat[3,i]=df
		nStartA=3+(0:6)*nCvOut
		jj=1:nCvOut
		confIntMat[nStartA[1]+jj,i]=Rmedian(cv.rsqList[[i]][[1]])
		confIntMat[nStartA[2]+jj,i]=Rmedian(cv.rsqList[[i]][[3]])
		confIntMat[nStartA[3]+jj,i]=Rmedian(FtoT(cv.rsqList[[i]][[3]]*df,1,df))
		confIntMat[nStartA[4]+jj,i]=Rmean(other.fstat < lowb05M)
 		confIntMat[nStartA[5]+jj,i]=Rmean(other.fstat < 0)
		confIntMat[nStartA[6]+jj,i]=Rmean(other.fstat[,ii05] < lowb05M[,ii05])
 		confIntMat[nStartA[7]+jj,i]=Rmean(other.fstat[,ii05] < 0)
 		}
 	colnames(confIntMat)=names(cv.rsqList)
 	confIntMat
 	}




##%%%%%%% Function to compute coverage probability for all cv,outL
## using means rather than medians

makeConfIntMatLCVM=function(cv.rsqList){
	nCvOut=nrow(cv.rsqList[[1]][[3]])
	CvOutNames=rownames(cv.rsqList[[1]][[3]])
	rsqNames=rownames(cv.rsqList[[1]][[1]])
 	confIntMat=matrix(0,3+(7*nCvOut),length(cv.rsqList))
      	 rownames(confIntMat)=c("mean in.Pr(F>0)","pr in.Pr(F)<0.05 =ii05",
			"degrees of freedom",
			paste("mean",rsqNames),
			paste("mean",CvOutNames),
			paste("mean tstat of",CvOutNames),
 			paste("pr",CvOutNames,"below 95% lowb"),
 			paste("pr",CvOutNames,"<0"),
 			paste("pr",CvOutNames,"below 95% lowb ii05"),
 			paste("pr",CvOutNames,"<0 ii05")
			)
 
 	for(i in 1:length(cv.rsqList)){
		medFstatProb=mean(cv.rsqList[[i]][[5]][6,])
 		ii05=cv.rsqList[[i]][[5]][6,]<=0.05
 		in.fstat=cv.rsqList[[i]][[5]][5,]
 		df=cv.rsqList[[i]][[5]]["df_modI",1]
 		other.fstat=cv.rsqList[[i]][[3]]*df
		lowb05=qf(0.05,1,df,in.fstat)
		lowb05[in.fstat<=lowb05]=0.0001
 		lowb05M=matrix(1,nCvOut,1)%*%lowb05
 
 		confIntMat[1,i]=medFstatProb
 		confIntMat[2,i]=sum(ii05)/length(ii05)
		confIntMat[3,i]=df
		nStartA=3+(0:6)*nCvOut
		jj=1:nCvOut
		confIntMat[nStartA[1]+jj,i]=Rmean(cv.rsqList[[i]][[1]])
		confIntMat[nStartA[2]+jj,i]=Rmean(cv.rsqList[[i]][[3]])
		confIntMat[nStartA[3]+jj,i]=Rmean(FtoT(cv.rsqList[[i]][[3]]*df,1,df))
		confIntMat[nStartA[4]+jj,i]=Rmean(other.fstat < lowb05M)
 		confIntMat[nStartA[5]+jj,i]=Rmean(other.fstat < 0)
		confIntMat[nStartA[6]+jj,i]=Rmean(other.fstat[,ii05] < lowb05M[,ii05])
 		confIntMat[nStartA[7]+jj,i]=Rmean(other.fstat[,ii05] < 0)
 		}
 	colnames(confIntMat)=names(cv.rsqList)
 	confIntMat
 	}


##%%%%%%% Function to collect all relvant data in an array for all cv,outL

makeConfIntArrLL=function(cv.rsqList){
	nCvOut=nrow(cv.rsqList[[1]][[3]])
	CvOutNames=rownames(cv.rsqList[[1]][[3]])
	rsqNames=rownames(cv.rsqList[[1]][[1]])
 	confIntArr=array(0,c(length(cv.rsqList),3+(5*nCvOut),200))
		nameS1=names(cv.rsqList)
      	nameS2=c("in.Pr(F>0)","in.Pr(F)<0.05",
			"degrees of freedom",rsqNames,CvOutNames,
			paste("tstat of",CvOutNames),
 			paste(CvOutNames,"below 95% lowb"),
 			paste(CvOutNames,"<0")
			)
		nameS3=colnames(cv.rsqList[[1]][[3]])
	dimnames(confIntArr)=list(nameS1,nameS2,nameS3)

 	for(i in 1:length(cv.rsqList)){
		FstatProb=(cv.rsqList[[i]][[5]][6,])
 		ii05=cv.rsqList[[i]][[5]][6,]<=0.05
 		in.fstat=cv.rsqList[[i]][[5]][5,]
 		df=cv.rsqList[[i]][[5]]["df_modI",]
 		other.fstat=cv.rsqList[[i]][[3]]*df
		lowb05=qf(0.05,1,df,in.fstat)
		lowb05[in.fstat<=lowb05]=0.0001
 		lowb05M=matrix(1,nCvOut,1)%*%lowb05
 
 		confIntArr[i,1,]=FstatProb
 		confIntArr[i,2,]=ii05
		confIntArr[i,3,]=df
		nStartA=3+(0:6)*nCvOut
		jj=1:nCvOut
		confIntArr[i,nStartA[1]+jj,]=cv.rsqList[[i]][[1]]
		confIntArr[i,nStartA[2]+jj,]=cv.rsqList[[i]][[3]]
		confIntArr[i,nStartA[3]+jj,]=FtoT(cv.rsqList[[i]][[3]]*df,1,df)
		confIntArr[i,nStartA[4]+jj,]=other.fstat < lowb05M
 		confIntArr[i,nStartA[5]+jj,]=other.fstat < 0
 		}
 
 	confIntArr
 	}






##%%%%%%%%%% Functions to compute r^2, and diff-r^2 tstats %%%%%%%%%%%

rsqVar=function(rsqStat,n,p){
	nuM=(4*rsqStat*((1-rsqStat)^2)*(n-p-1)^2)
	deN=(n^2-1)*(n+3)
	rsqvar=nuM/deN
	rsqvar
	}

rsqToT=function(rsqStat,n,p){
	rsqvar=rsqVar(rsqStat,n,p)
	rsqT=rsqStat/sqrt(rsqvar)
	rsqT
	}

rsqDtoT=function(rsqStat1,rsqStat2,n1,p1,n2,p2){
	rsqvar1=rsqVar(rsqStat1,n1,p1)
	rsqvar2=rsqVar(rsqStat2,n2,p2)
	rsqstdA=sqrt(rsqstd1+rsqstd2)
	rsqdelta=rsqStat1-rsqStat2
	rsqdeltaT=rsqdelta/rsqstdA
	rsqdeltaT
	}





##%%%%%%%%%%%%%%%%%%%%%% Bootstrap f^2, adj.f^2 %%%%%%%%%%%%%%%%%%%%%%%

makeFsqBootList=function(tar,inS,ImatX,ImatBoot,jjVec){
	jjList=(1:length(jjVec))
	names(jjList)=as.character(jjVec)
	fsqBootList=list()

    for(j in jjVec){
	Nmod=ncol(ImatX)
	Nboot=ncol(ImatBoot)
	jjj=1:ncol(inS)
	jj=setdiff(jjj,j)
	jjjcoef=1:(ncol(inS)+1)
	jjcoef=setdiff(jjjcoef,j+1)
	fsqMatBoot=matrix(0,ncol(ImatBoot),Nmod)
	adjFsqMatBoot=matrix(0,ncol(ImatBoot),Nmod)
	colnames(fsqMatBoot)=paste("mod",1:Nmod,sep="")
	colnames(fsqMatBoot)=paste("mod",1:Nmod,sep="")
	rsqMatBoot=fsqMatBoot
	adjRsqMatBoot=fsqMatBoot
	rsqMatBootD=fsqMatBoot
	adjRsqMatBootD=fsqMatBoot

	rownames(rsqMatBoot)=paste("r^2 sample",1:ncol(ImatBoot),sep="")
	rownames(adjRsqMatBoot)=paste("adj.r^2 sample",1:ncol(ImatBoot),sep="")
	rownames(rsqMatBootD)=paste("r^2D sample",1:ncol(ImatBoot),sep="")
	rownames(adjRsqMatBootD)=paste("adj.r^2D sample",1:ncol(ImatBoot),sep="")

	for(i in 1:Nmod){
		iiin=ImatX[,i];message(as.character(i))
			for(l in 1:Nboot){
				iib=iiin[ImatBoot[,l]]
				t1in=tar[iib]
				insI=inS[iib,]
				insD=inS[iib,jj]
				modD=lm(t1in ~ insD)
				modI=lm(t1in ~insI)
				modDsum=summary(modD)
				modIsum=summary(modI)
				rsqMatBoot[l,i]=modIsum$r.squared
				adjRsqMatBoot[l,i]=modIsum$adj.r.squared
				rsqMatBootD[l,i]=modDsum$r.squared
				adjRsqMatBootD[l,i]=modDsum$adj.r.squared	
				}
		}

	## Compute in and out of sample f^2
	fsqMatBoot=(rsqMatBoot-rsqMatBootD)/(1-rsqMatBoot) #NOTE!! Using (rsq_all) data
	adjFsqMatBoot=(adjRsqMatBoot-adjRsqMatBootD)/(1-adjRsqMatBoot) #NOTE!! Using (rsq_all) data
	rownames(fsqMatBoot)=paste("f^2 sample",1:ncol(ImatBoot),sep="")
	rownames(adjFsqMatBoot)=paste("adj.f^2 sample",1:ncol(ImatBoot),sep="")

	shortList=list()
	shortList[[1]]=rsqMatBoot
	shortList[[2]]=adjRsqMatBoot
	shortList[[3]]=rsqMatBootD
	shortList[[4]]=adjRsqMatBootD
	shortList[[5]]=fsqMatBoot
	shortList[[6]]=adjFsqMatBoot

	names(shortList)=c("rsqMatBoot","adjRsqMatBoot","rsqMatBootD",
				     "adjRsqMatBootD","fsqMatBoot","adjFsqMatBoot")
	J=jjList[as.character(j)]
	fsqBootList[[J]]=shortList
	}	
	fsqBootList
	}



##%%%%%%%%%%%% Formula based R^2 shrinkage estimators %%%%%%%%%%%%%%%

##%%%%%%%%%%%%%% Formula based estimators of outL.r.squaredI %%%%%%%%%%%%%
## In this section we define a number of formula based estimators that 
## are found in the review paper Raju et al 1997.

##%%%%%%%%%%%% Formulas for fixed X %%%%%%%%%%%%%%%%%
## Smith formula
rsqHatS=function(r.squared,N,p){
	ouT=1-(N/(N-p))*(1-r.squared)
	ouT}

## Wherry formula
rsqHatW=function(r.squared,N,p){
	ouT=1-((N-1)/(N-p))*(1-r.squared)
	ouT}

## Ezekiel formula
rsqHatE=function(r.squared,N,p){
	ouT=1-((N-1)/(N-p-1))*(1-r.squared)
	ouT}

## Olkin & Pratt formula
rsqHatOP=function(r.squared,N,p){
	term1=((N-3)*(1-r.squared))/(N-p-1)
	term2=1+(2*(1-r.squared))/(N-p+1)
	term3=(8*(1-r.squared)^2)/((N-p-1)*(N-p+3))
	ouT=1-term1*(term2+term3)
	ouT}

## Pratt approximation formula
rsqHatP=function(r.squared,N,p){
	term1=((N-3)*(1-r.squared))/(N-p-1)
	term2=1+(2*(1-r.squared))/(N-p-2.3)
	ouT=1-term1*term2
	ouT}

## Herzberg approximation formula
rsqHatH=function(r.squared,N,p){
	term1=((N-3)*(1-r.squared))/(N-p-1)
	term2=1+(2*(1-r.squared))/(N-p+1)
	ouT=1-term1*term2
	ouT}

## Claudy's empircally based formula
rsqHatC=function(r.squared,N,p){
	term1=((N-4)*(1-r.squared))/(N-p-1)
	term2=1+(2*(1-r.squared))/(N-p+1)
	ouT=1-term1*term2
	ouT}


##%%%%%%%% Formulas for random X %%%%%%%%%%%%%%%

## Lord and Nicholson formula
## Raju (11), Pin (17) called "Lord-2" there
rsqOutHatLN=function(r.squared,N,p){
	term1=(N+p+1)/N
	term2=((N-1)/(N-p-1))*(1-r.squared)
	ouT=1-term1*term2
	ouT}

## Darlington and Stein formula
## as given by Pin Thesis 1999 and Raju et al
## Raju (12), Pin (19)
rsqOutHatDS=function(r.squared,N,p){
	term1=(N+1)/N
	term2=(N-1)/(N-p-1)
	term3=((N-2)/(N-p-2))*(1-r.squared)
	ouT=1-term1*term2*term3
	ouT}


## Burket's formula (originally for correlation rather than R^2)
## Raju (13)^2 Pin (18)^2
rsqOutHatBU=function(r.squared,N,p){
	term1=(N*r.squared-p)^2
	term2=r.squared*(N-p)^2
	ouT=term1/term2
	ouT}


## Browne's formula with Wherry formula to unbias r.squared
## as found in Raju et al p.296-297
rsqOutHatBRW=function(r.squared,N,p){
	rsqUnBias=rsqHatW(r.squared,N,p)
	term1=((N-p-3)*rsqUnBias^2)+rsqUnBias
	term2=((N-2*p-2)*rsqUnBias)+p
	ouT=term1/term2
	ouT}


## Browne's formula with Olkin and Pratt formula to unbias r.squared
## as found in Raju et al p.296-297
rsqOutHatBROP=function(r.squared,N,p){
	rsqUnBias=rsqHatOP(r.squared,N,p)
	term1=((N-p-3)*rsqUnBias^2)+rsqUnBias
	term2=((N-2*p-2)*rsqUnBias)+p
	ouT=term1/term2
	ouT}

## We don't use the Claudy formulas found
## in Pin (22), (23) since they may be mistakes

## Rozeboom formula
rsqOutHatR=function(r.squared,N,p){
	ouT=1-((N+p)/(N-p))*(1-r.squared)
	ouT}



##%%%%%%%%%%%%% Function to compute all these at the same time %%%%%%%

rsqHatAll=function(r.sq,N,p){
	rsqH.S=rsqHatS(r.sq,N,p)
	rsqH.W=rsqHatW(r.sq,N,p)
	rsqH.E=rsqHatE(r.sq,N,p)
	rsqH.OP=rsqHatOP(r.sq,N,p)
	rsqH.P=rsqHatP(r.sq,N,p)
	rsqH.H=rsqHatH(r.sq,N,p)
	rsqH.C=rsqHatC(r.sq,N,p)
	rsqH.LN=rsqOutHatLN(r.sq,N,p)
	rsqH.DS=rsqOutHatDS(r.sq,N,p)
	rsqH.BU=rsqOutHatBU(r.sq,N,p)
	rsqH.BRW=rsqOutHatBRW(r.sq,N,p)
	rsqH.BROP=rsqOutHatBROP(r.sq,N,p)
	rsqH.R=rsqOutHatR(r.sq,N,p)
	nameS=paste("rsqH.",c("S","W","E","OP","P","H","C",
				  "LN","DS","BU","BRW","BROP","R"),sep="")
	ouT=matrix(0,length(r.sq),length(nameS))
	colnames(ouT)=nameS
	for(i in 1:length(nameS)){
		ouT[,i]=get(nameS[i])
		}
 	ouT
	}



##%%%%% Define a function to compute all the relevant rsq, cv.rsq, N, p etc. %%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
crossValA=function(taR,inS){
	Fit=lm(taR~inS)
	sumFit=summary(Fit)
	coeF=sumFit$coefficients
	Fit.p=sum(!is.na(coeF[2:nrow(coeF),"Estimate"]))
	Fit.N=nrow(inS)
	Fit.r.squared=sumFit$r.squared
	Fit.cv=crossValD1S(taR,inS)[1:3,]
	Fit.LN=rsqOutHatLN(Fit.r.squared,Fit.N,Fit.p)
	Fit.DS=rsqOutHatDS(Fit.r.squared,Fit.N,Fit.p)
	Fit.R=rsqOutHatR(Fit.r.squared,Fit.N,Fit.p)
	## Collect all of these
	ouT=matrix(c(Fit.N,Fit.p,Fit.cv,Fit.LN,Fit.DS,Fit.R),8,1)
	ouT[1:2,]=round(ouT[1:2,])
	colnames(ouT)="Fit"
	rownames(ouT)=c("N","p","r.squared","adj.r.squared","cv.r.squared", 
				       "LN.r.squared","DS.r.squared","R.r.squared")
	ouT}


## Function to compute all but cv.r.squared using just r.squared, N, p
crossValAnoCV=function(rsq,N,p){
	Fit.p=p
	Fit.N=N
	Fit.adj.r.squared=adjRsq(rsq,N,p)
	Fit.cv=NA
	Fit.LN=rsqOutHatLN(rsq,Fit.N,Fit.p)
	Fit.DS=rsqOutHatDS(rsq,Fit.N,Fit.p)
	Fit.R=rsqOutHatR(rsq,Fit.N,Fit.p)
	ouT=matrix(c(Fit.N,Fit.p,rsq,Fit.adj.r.squared,NA,Fit.LN,Fit.DS,Fit.R),8,1)
	ouT[1:2,]=round(ouT[1:2,])
	colnames(ouT)="Fit"
	rownames(ouT)=c("N","p","r.squared","adj.r.squared","cv.r.squared", 
				       "LN.r.squared","DS.r.squared","R.r.squared")
	ouT}




##%%%%%%%%%%%%% Functions for variable level effect size %%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%% Formula based estimators of outL.r.squaredI %%%%%%%%%%%%%
## In this section we define a number of formula based estimators that 
## are found in the review paper Raju et al 1997.

##%%%%%%%%%%%% Formulas for fixed X %%%%%%%%%%%%%%%%%
## Smith formula
rsqHatS=function(r.squared,N,p){
	ouT=1-(N/(N-p))*(1-r.squared)
	ouT}

## Wherry formula
rsqHatW=function(r.squared,N,p){
	ouT=1-((N-1)/(N-p))*(1-r.squared)
	ouT}

## Ezekiel formula
rsqHatE=function(r.squared,N,p){
	ouT=1-((N-1)/(N-p-1))*(1-r.squared)
	ouT}

## Olkin & Pratt formula
rsqHatOP=function(r.squared,N,p){
	term1=((N-3)*(1-r.squared))/(N-p-1)
	term2=1+(2*(1-r.squared))/(N-p+1)
	term3=(8*(1-r.squared)^2)/((N-p-1)*(N-p+3))
	ouT=1-term1*(term2+term3)
	ouT}

## Pratt approximation formula
rsqHatP=function(r.squared,N,p){
	term1=((N-3)*(1-r.squared))/(N-p-1)
	term2=1+(2*(1-r.squared))/(N-p-2.3)
	ouT=1-term1*term2
	ouT}

## Herzberg approximation formula
rsqHatH=function(r.squared,N,p){
	term1=((N-3)*(1-r.squared))/(N-p-1)
	term2=1+(2*(1-r.squared))/(N-p+1)
	ouT=1-term1*term2
	ouT}

## Claudy's empircally based formula
rsqHatC=function(r.squared,N,p){
	term1=((N-4)*(1-r.squared))/(N-p-1)
	term2=1+(2*(1-r.squared))/(N-p+1)
	ouT=1-term1*term2
	ouT}


##%%%%%%%% Formulas for random X %%%%%%%%%%%%%%%

## Lord and Nicholson formula
## Raju (11), Pin (17) called "Lord-2" there
rsqOutHatLN=function(r.squared,N,p){
	term1=(N+p+1)/N
	term2=((N-1)/(N-p-1))*(1-r.squared)
	ouT=1-term1*term2
	ouT}

## Darlington and Stein formula
## as given by Pin Thesis 1999 and Raju et al
## Raju (12), Pin (19)
rsqOutHatDS=function(r.squared,N,p){
	term1=(N+1)/N
	term2=(N-1)/(N-p-1)
	term3=((N-2)/(N-p-2))*(1-r.squared)
	ouT=1-term1*term2*term3
	ouT}


## Burket's formula (originally for correlation rather than R^2)
## Raju (13)^2 Pin (18)^2
rsqOutHatBU=function(r.squared,N,p){
	term1=(N*r.squared-p)^2
	term2=r.squared*(N-p)^2
	ouT=term1/term2
	ouT}


## Browne's formula with Wherry formula to unbias r.squared
## as found in Raju et al p.296-297
rsqOutHatBRW=function(r.squared,N,p){
	rsqUnBias=rsqHatW(r.squared,N,p)
	term1=((N-p-3)*rsqUnBias^2)+rsqUnBias
	term2=((N-2*p-2)*rsqUnBias)+p
	ouT=term1/term2
	ouT}


## Browne's formula with Olkin and Pratt formula to unbias r.squared
## as found in Raju et al p.296-297
rsqOutHatBROP=function(r.squared,N,p){
	rsqUnBias=rsqHatOP(r.squared,N,p)
	term1=((N-p-3)*rsqUnBias^2)+rsqUnBias
	term2=((N-2*p-2)*rsqUnBias)+p
	ouT=term1/term2
	ouT}

## We don't use the Claudy formulas found
## in Pin (22), (23) since they may be mistakes

## Rozeboom formula
rsqOutHatR=function(r.squared,N,p){
	ouT=1-((N+p)/(N-p))*(1-r.squared)
	ouT}




##%%%%%%%%%%%%% Function to compute all these at the same time %%%%%%%

rsqHatAll=function(r.sq,N,p){
	rsqH.S=rsqHatS(r.sq,N,p)
	rsqH.W=rsqHatW(r.sq,N,p)
	rsqH.E=rsqHatE(r.sq,N,p)
	rsqH.OP=rsqHatOP(r.sq,N,p)
	rsqH.P=rsqHatP(r.sq,N,p)
	rsqH.H=rsqHatH(r.sq,N,p)
	rsqH.C=rsqHatC(r.sq,N,p)
	rsqH.LN=rsqOutHatLN(r.sq,N,p)
	rsqH.DS=rsqOutHatDS(r.sq,N,p)
	rsqH.BU=rsqOutHatBU(r.sq,N,p)
	rsqH.BRW=rsqOutHatBRW(r.sq,N,p)
	rsqH.BROP=rsqOutHatBROP(r.sq,N,p)
	rsqH.R=rsqOutHatR(r.sq,N,p)
	nameS=paste("rsqH.",c("S","W","E","OP","P","H","C",
				  "LN","DS","BU","BRW","BROP","R"),sep="")
	ouT=matrix(0,length(r.sq),length(nameS))
	colnames(ouT)=nameS
	for(i in 1:length(nameS)){
		ouT[,i]=get(nameS[i])
		}
 	ouT
	}


adjRsq=function(rsq,N,p){
	ouT=1-((1-rsq)*(N-1))/(N-p-1)
	ouT}

##%%%%% Define a function to compute all the relevant rsq, cv.rsq, N, p etc. %%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
crossValA=function(taR,inS){
	Fit=lm(taR~inS)
	sumFit=summary(Fit)
	coeF=sumFit$coefficients
	Fit.p=sum(!is.na(coeF[2:nrow(coeF),"Estimate"]))
	Fit.N=nrow(inS)
	Fit.r.squared=sumFit$r.squared
	Fit.cv=crossValD1S(taR,inS)[1:3,]
	Fit.LN=rsqOutHatLN(Fit.r.squared,Fit.N,Fit.p)
	Fit.DS=rsqOutHatDS(Fit.r.squared,Fit.N,Fit.p)
	Fit.R=rsqOutHatR(Fit.r.squared,Fit.N,Fit.p)
	## Collect all of these
	ouT=matrix(c(Fit.N,Fit.p,Fit.cv,Fit.LN,Fit.DS,Fit.R),8,1)
	ouT[1:2,]=round(ouT[1:2,])
	colnames(ouT)="Fit"
	rownames(ouT)=c("N","p","r.squared","adj.r.squared","cv.r.squared", 
				       "LN.r.squared","DS.r.squared","R.r.squared")
	ouT}

## Function to compute all but cv.r.squared using just r.squared, N, p
crossValAnoCV=function(rsq,N,p){
	Fit.p=p
	Fit.N=N
	Fit.adj.r.squared=adjRsq(rsq,N,p)
	Fit.cv=NA
	Fit.LN=rsqOutHatLN(rsq,Fit.N,Fit.p)
	Fit.DS=rsqOutHatDS(rsq,Fit.N,Fit.p)
	Fit.R=rsqOutHatR(rsq,Fit.N,Fit.p)
	ouT=matrix(c(Fit.N,Fit.p,rsq,Fit.adj.r.squared,NA,Fit.LN,Fit.DS,Fit.R),8,1)
	ouT[1:2,]=round(ouT[1:2,])
	colnames(ouT)="Fit"
	rownames(ouT)=c("N","p","r.squared","adj.r.squared","cv.r.squared", 
				       "LN.r.squared","DS.r.squared","R.r.squared")
	ouT}


## Define minShrTstat, avgShrTstat
minShrTstat=function(shrVec){
	shrN=max(shrVec[!is.na(shrVec)])
	ouT=0
	if(shrN>0) ouT=sqrt(shrN)
	ouT}

avgShrTstat=function(shrVec){
	shrN=mean(shrVec[!is.na(shrVec)])
	ouT=0
	if(shrN>0) ouT=sqrt(shrN)
	ouT}


